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Abstract. Inspired by the physics of magnetohydrodynamics (MHD) a simplified 
coupled Burgers- like model in one dimension (Id), a generalization of the Burgers 
model to coupled degrees of freedom, is proposed to describe lcflVIHD. In addition 
to MHD, this model serves as a Id reduced model for driven binary fluid mixtures. 
Here we have performed a comprehensive study of the universal properties of the 
generalized ci-dimensional version of the reduced model. We employ both analytical 
and numerical approaches. In particular, we determine the scaling exponents and 
the amplitude-ratios of the relevant two-point time-dependent correlation functions 
in the model. We demonstrate that these quantities vary continuously with the 
amplitude of the noise cross-correlation. Further our numerical studies corroborate 
the continuous dependence of long wavelength and long time-scale physics of the 
model on the amplitude of the noise cross-correlations, as found in our analytical 
studies. We construct and simulate lattice-gas models of coupled degrees of freedom 
in Id, belonging to the universality class of our coupled Burgers-like model, which 
display similar behavior. We use a variety of numerical (Monte-Carlo and Pseudo- 
spectral methods) and analytical (Dynamic Renormalization Group, Self-Consistent 
Mode-Coupling Theory and Functional Renormalization Group) approaches for our 
work. The results from our different approaches complement one another. Possible 
realizations of our results in various noncquilibrium models are discussed. 
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1. Introduction 

Physical descriptions of many natural driven systems involve coupled dynamics of several 
degrees of freedom. A prominent example is the dynamics of a driven symmetric mixture 
of a miscible binary fluid pQ. Coarse-grained dynamical descriptions of such a system are 
in terms of the local velocity field v(r, t) and the difference in the local concentrations of 
the two components ip(r,t). The feature which sets this system apart from a passively 
advected system is that density fluctuations, while being advected by the flow, create 
stresses which in turn feed back on the flow. Thus, this forms an example of active 
advection. Typical experimental measurements include correlation functions of the 
dynamical variables and their universal scaling properties [2] in the non-equilibrium 
statistical steady states (NESS). A similar system with a closely related theoretical 
structure is magnetohydrodynamics in three dimensions (3dMHD). It deals with the 
coupled evolution of the velocity v(r, t) and magnetic b(r, t) fields in quasi-neutral 
plasmas [3j HI [5]. There the v-field advects and stretches the b-field lines; the b- 
field in turn acts on the v-field through the Lorentz force. As in driven binary fluids, 
one is typically interested in the correlation functions and their universal properties 
characterized by appropriate scaling exponents in the NESS. Despite the differences in 
the microscopic physics of these two systems, they have a great many commonalities in 
the physical descriptions at large scales and long times. In particular, they show similar 
universal scaling behavior in the hydrodynamic limit. Often one asks similar questions 
concerning scaling behavior in the NESS and the nature of the dynamical equations these 
systems follow. Inspired by these generalities a Id generalized Burgers model [6] was 
proposed as a simpler reduced model, which captures many qualitative features of two 
real systems we described above. This model consists of the Burgers equation for the 
velocity field, containing a coupling term to the second field (representing feedback), 
which in turn obeys a diffusion-advection equation. We generalize this model to d- 
dimensions, where it has the same continuous symmetries as the binary fluid mixture 
model or 3<iMHD model equations, and investigate its universal properties as a function 
of the dimension d. 

In the vicinity of a critical point, equilibrium systems show universal scaling 
properties for thermodynamic functions and correlation functions. These are 
characterized by universal scaling exponents that depend on the spatial dimension 
d and the symmetry of the order parameter (e.g., Ising, XY etc.) [7j, but not on 
the material parameters that specify the (bare) Hamiltonian. Prominent exceptions 
are the 2d XY model and its relatives, where the renormalization group flow is 
characterized by a fixed line and as a consequence the scaling exponents exhibit a 
continuous dependence on the stiffness parameter. Equilibrium dynamics close to critical 
points also show universality though the dynamic scaling exponents, which characterize 
the time-dependence of unequal time correlation functions, now also depend upon the 
presence or absence of conservation laws and the non-dissipative terms in the dynamical 
equations [8]. A different situation arises in driven, dissipative, nonequilibrium systems 
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with NESS. Significant advances have been made in classifying the physics of non- 
equilibrium systems at long time and large length scales into universality classes. It has 
been shown that standard universality classes in critical dynamics are quite robust to 
detailed-balance violating perturbations [9]. Novel features are found only for models 
with conserved order parameter and spatially anisotropic noise correlations. In contrast, 
truly non-equilibrium dynamic phenomena, whose steady states cannot be described 
in terms of Gibbsian distributions, are found to be rather sensitive to all kinds of 
perturbations. Prominent examples are driven diffusive systems [TO], and fluid- and 
magnetohydrodynamic- turbulence [U [HI H2]- In contrast to systems in equilibrium, 
what characterizes the universality classes in nonequilibrium systems remains a yet 
unsolved question. 

In this article we examine the particular issue of universality in the non-equilibrium 
in the context of driven coupled generalized Burgers model [6] in ^-dimensions. We 
show that its NESS depends sensitively on the parameters of the model. A brief 
account of the physics of this model has been reported in Ref.[13]. Here we extend 
the results of Ref . [13] , and present a comprehensive study of the universal properties 
of the model. In order to study our model systematically, we have employed a variety 
of analytical and numerical techniques, all of which together bring out a coherent and 
consistent picture: We find that the universal properties of the model characterized by 
dimensionless amplitude ratios and scaling exponents of various correlation functions 
in the model depend explicitly on the noise crosscorrelations. Our results provide 
valuable insight on the issue of universality in coupled nonequilibrium systems and, 
in particular, what might characterize universality classes in a simple coupled model 
discussed here. We use analytical methods including dynamic renormalization group 
(DRG), self-consistent mode coupling methods (SCMC) and functional renormalization 
group (FRG) to calculate the relevant scaling exponents and amplitude-ratios of the 
correlation functions. Furthermore, to complement our analytical results we use pseudo- 
spectral methods to numerically solve the stochastically driven model equations in one 
(Id) and two (2d) space dimensions and use them to calculate scaling exponents and 
amplitude ratios. We, in addition, construct lattice-gas models in one space dimension 
belonging to the same universality class as the model equations of [13], following 
closely the lattice-gas models for the Kardar-Parisi- Zhang (KPZ) equation for growing 
surfaces [H]. We perform Monte-Carlo simulations on this coupled lattice-gas model and 
calculate the scaling exponents and the amplitude ratios. Our analytical and numerical 
results agree well with each other. The organization of the rest of the paper is as follows: 
In Section 12.11 we consider the model equations and the noise correlations which we use 
in our work here. In Sec. 12.21 we describe the noise correlations which we use in our 
analytical work. Then in the Section [3] we describe our constructions of the lattice- 
gas models in Id belonging to the universality class of the continuum model equations. 
Then in Sections H] and [5] we present our analytical and numerical results respectively, 
concerning the scaling exponents and the amplitude-ratios. We finally summarize our 
work in Secj6j 
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2. The stochastically driven model equations 

2.1. Construction of the model equations 

Let us begin by recalling the general principles which go in setting up the one- 
dimensional (Id) Burgers [15] as a reduced model for the 3d Navier-Stokes (NS) equation 
for the velocity field v. The NS equation is given by 
d~v 

— + v ■ Vv = -Vp/p + v v V\ + f„, (1) 

where p is the pressure, p is the density, v v is the kinematic viscosity and f v is an external 
force (deterministic or random). Equation (I5.2.2P has two constants of motion in the 
inviscid, unforced limit (y v = 0, i v = 0): (i) The kinetic energy J d d xv 2 and (ii) helicity 
J Av ■ V x v. Further equation ([1]) is invariant under the Galilean transformation 
and, as a consequence, is of the conservation law form. The Id simplest non-linear 
equation which obeys the conservation of kinetic energy (helicity cannot be defined in 
Id properly, since it involves curl of a vector) and also is Galilean invariant is the famous 
Burgers equation [TS] 

dv 

757 + ^dxv 2 = vd xx v + /. (2) 
at 2 

Here v is the Id Burgers velocity field - it is a pressureless fluid since the pressure has 
been dropped, v is the Burgers viscosity and / is an external force. 

Fluid systems having coupled degrees of freedom require descriptions more general 
than pure one-component fluid systems. Two such well-known examples are 



(i) Symmetric (50-50) binary fluid mixture:- The coarse-grained dynamics of a driven 
symmetric binary fluid mixture is described by two continuum variables: a velocity 
field v(r, t) and a concentration-gradient field b(r, t) = V^>(r, i) where ip is the relative 
difference of local densities of the two components [1]. The coupled equations of motion, 
in the incompressible limit of the dynamics, consist of the generalized Navier-Stokes 
(which includes the stresses from the concentration gradients) for v and a diffusion- 
advection equation for b 

v + v-Vv = a bV ■ b + z/V 2 v + f v , (3) 

at p 

^_ + V(vb)= /U V 2 b + f b , (4) 

together with the incompressibility condition V • v = 0. We also have a > for 
thermodynamic stability; p, is diffusivity of the concentration gradient. Functions f v and 
f{> are stochastic forces required to maintain a statistical steady-state. Force f& is curl- 
free. Equations ([3]) and (j4j) are Galilean invariant and admit the following conservation 
laws in the inviscid, unforced limit [I]: 

• Spatial integral of the square of the concentration ip: C = ~ J d 3 r[ip(r, t} 2 , 

• Total energy E = ~ j d 3 r[v(r,t)\ 2 + ab(r, t) 2 }. 
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We now ask for a Id model, whose relation with equations ([3]) and (jlj) is same as that 
between the Burgers equation and the Navier-Stokes equation. The simplest of such Id 
equations, to the leading orders in the gradients and bilinear order in the fields u and 
b are [6] 

^ + ^d x u 2 + ^d x b 2 = vd xx u + /, (5) 
db 

— + d x (ub) =fxd xx b + g. (6) 

Here, u and b are the Id Burgers velocity and concentration gradient fields, v and p are 
diffusion coefficients for u and b, and / and g are external sources. Equations (jSJ) and 
(EJ) the Id analogs of E and C as defined above in the inviscid {y — p — 0), unforced 
(/ = g — 0) limit obey the two conservation laws mentioned above. Note that equations 
([5]) and ([6]) are invariant under the Galilean transformation. 

Magnetohydrodynamics in three dimensions:- The subject of Magnetohydrodynamics 
(MHD) [3], HI [5] deals with the dynamics of a quasi- neutral plasma in terms of coupled 
evolution of velocity v(r, t) and magnetic b(r, t) fields, the equations of motion consist 
of the generalized Navier-Stokes equation containing the stresses from the magnetic 
fields for the v field and the Induction equation for the b field: 

<9v „ Vp (V x b) x b - 

+ v . V v = + ± -± + vV\ + f v , 7 

at p Anp 

db 

- + Vx(vxb)= /1 V 2 b + f, (8) 

These are to be supplemented by V ■ b = 0, and, in case of an incompressible fluid, 
V ■ v = 0. Here, u, p are kinetic and magnetic viscosities and i v , fb are the external 
sources, with fj, being solenoidal. Equations (J7|) and (jHJ) are Galilean invariant and 
admit the following conserved quantities in the inviscid and unforced limit: (i) Total 
energy (the sum of the kinetic and magnetic energies) E = | J d 3 r[v 2 + b 2 }, and (ii) 
Crosshelicity C = J d 3 r v • b As in the previous example of binary fluid mixture we 
what the Id model should be whose relation with the 3<iMHD equation above is same 
as that of Burgers with Navier-Stokes. Such Id equations are the same as (jSJ) and (jSJ), 
which conserves the Id analog of the total crosshelicity as well Further, Eqs. of motion 
and (jHJ) are invariant under the Galilean transformation. Note also that under parity 
reversal r — > — r, v and b transform differently: v — » — v and b — > b. 

The rf-dimensional generalization [13] of the ld-model equations (jSJ) and ([6]) are, 



^ + ^Vu 2 + ^Vb 2 = , V 2 u + f, (9) 



db 

~dt 



A 3 V(u-b) = /i V 2 b + g. (10) 



| The third invariant of 3dMHD, i.e., the total magnetic helicity cannot be properly defined in Id, 
since it involves curl of a vector. Hence it is not considered in writing down the Id model. 
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We refer to these equations as the Generalized Burgers Model (henceforth GBM). 
Here, u and b are <i-dimensional Burgers fields, respectively. Parameters Ai, A2 and 
A3 are actually identical to unity but kept for formal book keeping purposes in the 
renormalization group and mode coupling approaches below. Parameters u and ji^ 
are bare (unrenormalized; see below) viscosities. Functions f and g are external noise 
sources required to maintain a statistical steady state. The noise g, in general, may 
have both a non-zero divergence and a curl (see the discussion just before Eq. (|45p . The 
quantities of interests are the correlation functions as functions of Fourier wavevector k 
and frequency uj 

C%{k,U>) = (Ui(k,u)Uj(-k,-U>)), (11) 

(7*(k,u;) = (6 i (k,a;)6 i (-k,-u;)>, (12) 
C^ j {k i u) = {u i (k 1 u)b j (-k i -u)). (13) 

From the properties of the fields u and b under the reversal of parity as discussed above 
Cf, and C. 6 ,- are even functions of k and Cf, is an odd function of k. 

Upon introducing a new set of fields, u = Vh and b = V0, the GBM [6] maps onto 
a model for drifting lines introduced by Ertas and Kardar [16] 



f + y (V/*) 2 + y (V0) 2 = u V 2 h + B U (14) 

^ + X 3 (Vh) • (V0) = /i o V 2 + 9 2 , (15) 

where fields h and <fi are the local, instantaneous longitudinal and transverse fluctuations, 
respectively, around the mean position of the drifting line. The functions 9\ and 82 are 
given by f = V#i and g = V9 2 . If the fluid and magnetic viscosities are equal, v$ = /io, 
the Burgers-like MHD model, Eqs. (Q and ffTUl) . can also be mapped onto a model of 
coupled growing surfaces. Further by introducing Elsasser variables, z^ = u ± b one 
finds (setting Ai = A2 = A3 = 1), for 770 = is = fiQ, 

^ + \^z± 2 = r /0 V 2 z ± + 0± , (16) 

where Of = (f t ± g t )/2. 

Each of the Elsasser variables obeys a stochastically driven Burgers equation, where 
the coupling between the two fields arises solely due to the cross-correlations in 
the external (stochastic) driving forces. Physically, this model describes the growth of 
surfaces on two interpenetrating sub-lattices (say, A and B), where the growth of each of 
the surfaces follows a KPZ dynamics. The dynamics becomes coupled since depositions 
of particles on the A and B sub-lattices are correlated with each other. Such a coupled 
surface growth problem can be mapped to a related equilibrium problem of a pair of 
two Directed Polymers (DP) in random medium. A DP in d + 1 dimension is just 
a directed string stretched along one particular direction with free fluctuations in all 
other d transverse directions. When placed in a random medium competitions between 
the elastic energy of the string and the random potential of the medium lead to phase 
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transitions between smooth and rough phases [17] . The DP phases and phase transition 
can be mapped exactly to the phases and phase transitions in the KPZ surface growth 
model. In the present case, the pair of variables model the free energies of the two 
DPs in a random medium. The variable x refers to the directions transverse to the 
DP and t is the longitudinal dimension. However, the external noise sources 6^ are in 
general cross-correlated leading to interesting phase diagram of the coupled system [18] 
(see also below). 



2.2. Noise distributions 

The noise sources f and g or alternatively 6 ± are chosen to be Gaussian-distributed with 
zero mean and specified variances. It should be noted that, in nonequilibrium situations 
there are no restrictions on the noise variances as there are no detailed balance conditions 
relating the diffusivities and the noise- variances, unlike for systems in equilibrium |19j . 
Furthermore, for analytical conveniences we assume them to be conserved noise sources 
which suffices for our purposes in this article. In particular we make the following 
choices: 

(/i(k,t)/ i (-k,0)> = 2k l k j D^(k)6(t), (17) 
(^(k,%(-k,0)> = 2k i k j Di°\k)6(t), (18) 
(/ l (k,%(-k,0)) = 2z^. (0) (k)5(t), (19) 

where the functions D^(k) are even and D^ \k) is odd in k, respectively. A 
superscript refers to the bare (unrenormalized) quantities. The above structures 
of the auto-correlators of /, and g\ are the simplest choices consistent with their 
tensorial structures and the conservation law form of the equations of motion (Q and 
f fTOl) . Equations ( JT7l) and (TIB"]) are invariant under inversion, rotation and exchange 
of % with j. We take the noise cross-correlation, equation ( |T9l) to be invariant under 
spatial inversion, but we allow it to break (i) rotational invariance, and (separately) 
(ii) symmetry with respect to an interchange of the cartesian indices i and j. The 
matrix D^°\\l) in general has a symmetric part Df- (k) s and an antisymmetric part 
Df^°\]s.) a with respect to interchanges between i and j. In this article most of our 

results involve a finite Df-^Ck) s but zero ££ (0) (k) o , although some effects of Df- 
are also discussed. Effects of the noise cross-correlations of the forms given in Eq. (jT9l) 
on the universal properties of the model Eqs. ([9]) and ( TTUT) are discussed briefly in 
Ref . [13] . The properties of noise correlators described above use the explicit symmetries 
of the GBM model and the fields u and b. Such choices, however, leave the functional 
forms (as functions of k) of Du\k), D^\\s.) and Df- (k), which are the amplitudes 
of the noise variances, arbitrary. In order to define the model completely by making 
specific choices, one further defines amplitudes of the symmetric and antisymmetric 
parts of the cross- correlations by the relations Z^. (0) (k) s = 2L>°(k)A; i £; i , |D°(k)| = D° s 
and Df^°\k) a Df^°\k) a = AD'f'k 4 ' where D° s and D® are amplitudes of the respective 
parts. Further, we choose Du\k) = D® and D^\k) = D®. 
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3. Coupled lattice-gas models in one dimension 

Studies of lattice-gas models to understand the long-time, large-scale properties of 
continuum model equations in non-equilibrium statistical mechanics have a long history. 
For example, several lattice-gas model for the KPZ equation have been proposed and 
studied (see, e.g., [201 El])- Often numerical simulations are used to investigate the 
scaling properties. Such studies have several advantages. The main advantage is: Due to 
the analytically intractable nature of such models it is often much easier to numerically 
simulate a lattice-gas model than to obtain the numerical solution of the corresponding 
noise-driven continuum model equation. In addition such studies allow us to explore and 
understand the applicability of the concept of universality classes, originally developed 
in the context of critical phenomena and equilibrium critical dynamics [7], in physical 
situations out of equilibrium, by comparing the lattice model results with the results 
from the continuum model. 

In this section we propose a lattice-gas model for the one-dimensional (Id) GBM, 
Eqs. ([9]) and (flOl . Our starting point is the observation made in Ref. [13], that in the 
hydrodynamic limit the effective (renormalized) Prandtl number P m = u/fi = 1 for the 
model Eqs. (J9J) and (fTUl) |§| Henceforth, since we are interested in the asymptotic scaling 
properties, we may set the bare magnetic Prandtl number = [Iq/vq = 1, i.e., take 
the two bare viscosities as identical. Upon setting our book keeping parameters to their 
physical values, Ai = A2 = A3 = 1, and introducing the Elsasser variables z = u ± b, 
one obtains a set of coupled Burgers equations 

d t z± + \d x (z±) 2 = uoftz* + e± , (20) 

where the coupling between z + and z~ is mediated only via the cross-correlations in the 
noise 8± = |(/ ± g) only. With the standard mapping, z + = d x h\, z~ = d x fi2, we may 
rewrite Eq. (1201) as a set of coupled KPZ equations for the height fields h± and h>2 

dth + \^{d x hi) 2 = u d xx hi + tpi, 

d t h 2 + ^{d x h 2 ) 2 = v d xx h 2 + tp 2 , (21) 

where the noises d x ipi t 2 = 0±. Further the auto-correlations of the fields are simply 
related to those of h and <fi: hi 2 — h± <fi. We write them in the real space 

(h(x,t)h(x',t')) = ^[(h 1 (x,t)h 1 (x',t')) + (h 2 (x,t)h 2 (x',t')) 

+ (hx(x, t)h 2 (x, t)) + (hix', t')h 2 (x, t)}, (22) 

(<f){x,t)<f){x',t')) = -[(h l (x,t)h l (x',1?)) + (h 2 (x,t)h 2 (x',t')) 

- (h x {x,t)h 2 {x' ,t')) - (h x (x' ,t')h 2 (x,t))\. (23) 

Constructions of lattice models for a single KPZ equation is well-documented in 
the literature [20l [21]. In such models the underlying space (substrate) is taken to be 

§ Technically, as discussed in Ref. [13J . renormalized P m = 1 is the fixed point of the model. 
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discrete. Particles are deposited randomly over the discrete lattice and they settle on 
different lattice points according to certain model-dependent local growth rules [20], [2Tj . 
In our case since each of the Eqs. (1211) has the structure of a KPZ equation, such growth 
rules can be used to represent each of Eqs. ([2~T]) . representing two growing surfaces over 
two sub-lattices. In such lattice models stochasticity enters into the model through 
the randomness in the deposition process. In our case we model noise-crosscorrelations 
[see Eqs. (I24H26I) below] by cross-correlating the randomness in depositions in the two 
sublattices. In the Fourier space the correlations of the noise sources ipi? have the form 



Note that though variances in Eqs. (I24H25P in real space are proportional to 8{x— x'), the 
third one is proportional to l/(x — x') (see Appendix 19.21) . Such noises can be generated 
by first generating two independent, zero- mean short range, white in space and time 
Gaussian noises, and then taking appropriate linear combinations of them. A numerical 
scheme to generate noises with variances ( [2411261) is given in Appendix ( 19. 4ft . Eqs. ( ]2TT) 
may be viewed as a coupled growth model describing the dynamics of two growing 
surfaces of two different particles on two interpenetrating sub-lattices A and B where 
the deposition process (represented by the noise sources ipi^) of the particles A and B 
are correlated but the local dynamics of the particles on the sublattices are independent 
of each other. In particular each sub-lattice has a local KPZ dynamics. Monte-Carlo 
simulations of our coupled lattice models yields fields hi and h 2 . The dynamics of the 
two growing surfaces on the two sublattices are coupled through the noise sources. The 
auto-correlators of h and <fi are then calculated by using relations ([23]) . 

We implement the Newman-Bray (hereafter NB) [22] and the Restricted Solid on 
Solid (RSOS) [20] algorithms for each of the KPZ equations. In the RSOS algorithm 
[20] noise sources with correlations characterized by Eqs. (I24H26I) are used in the 
random determination of the lattice sites which are to be updated in a given time- 
step; see below for more details. Since tpi^ are cross-correlated the sites of sub-lattices 
A and B which are updated in a given time-step get cross-correlated, which in turn 
induces cross-correlations in the height fields hi and h 2 . We have performed Monte- 
Carlo simulations on our coupled lattice model based on the RSOS lattice model for 
the KPZ dynamics. The results from the Monte-Carlo simulations of our coupled RSOS 
lattice model are discussed in Section [5. 2. II The details of generation of random numbers 
obeying correlations ( ]26l) are discussed in the Appendix ( 19.41 ). 

In the NB lattice-gas model, the mapping between the KPZ surface growth model 
and the equilibrium problem of directed polymer in a (quenched) random medium is 
exploited [TT] [22]. We extend this idea in our construction of a lattice-gas model for the 
Eqs. ( 1211 based on the NB algorithm which is equivalent to considering two directed 
polymers in a random medium. The two free-energies of the two polymers would then 
represent the heights of the two interpenetrating sub-lattices as described above. In this 



M)> 
M)> 



2D 5(t), 
2iD k/\k\5(t). 



(24) 
(25) 
(26) 
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coupled lattice-gas model, the noise cross-correlations represent effective interactions 
induced by the randomness of the embedding medium between the two polymers. The 
detailed numerical results from the model are presented in Section 15.2.21 

4. Field theory analysis 

Our field theoretic analytical studies include one-loop dynamic renormalization group 
(DRG), one-loop self-consistent mode coupling (SCMC) and functional renormalization 
group (FRG) studies. From previous studies on the closely related KPZ model for 
surface growth one has learned that DRG schemes are well-suited to study the scaling 
properties of the rough phase in Id and the smooth-to-rough phase transition in spatial 
dimensions d = 2 + e, e>0 [23]. In contrast, the SCMC and the FRG approches are 
known to yield results on the scaling properties of the rough phases in d = 1 dimensions 
and higher pH E3]. In d = 1 the results from the DRG and the SCMC/FRG schemes 
are identical |25j . 

We are interested in the physics in the scaling limit, i.e., at long time and length 
scales. In that limit the time-dependent two-point correlation functions are written in 
terms of the dynamic exponent z and the two roughness exponents Xh an d X(j> a ^jjl 

CW(x,t) = (/i(x,t)/i(0,0)) = x 2 ^f h (t/x z ) , (27) 
CW(x, t) = (0(x, t)0(O, 0)) = x 2 ^U(t/x z ) , (28) 
C x (x, t) = (h(pc, t)0(O, 0)) = sgn(x)x^ + ^/ x (t/x z ) . (29) 

Here, angular brackets () means averaging over the noise distributions. Functions fh, f<f, 
and f x are scaling functions of the scaling variable t/x z . Since C x (x) is an odd function 
of x, a signum function appears, sgn(x) = — sgn(— x). Ward identities resulting from the 
Galilean invariance of the model Eqs. (Q and fflUj) imply that Xh an d X<t> are identical 
[see below; see also Appendix 19. lj . Henceforth, we write Xu — Xh~ 1 — X<t>~ 1 = Xb = X- 
Clearly then the ratios of the various equal-time correlators are dimensionless numbers. 
One also defines widths 

W h (t) = V([Mx,t)-(fr(x,t))] 2 ) = y/ [(Mx,t) 2 )-(M"x^)F] (30) 

w&) = vfM 7 **) = VMx,m - (tfrt))*] , (3i) 

These are related to the two-point correlators measured at the same space and time. 
They exhibit growing parts for small t and yield the ratio Xh/ Z '- Wh(t), W^it) ~ t Xh ' z , 
and saturated parts for large t yielding the exponent Xh- Wh{t), W^{t) ~ L Xh for large 
t where L is the system size. The ratios of the amplitudes of the correlation functions 
or the widths in the steady-state yield the amplitude-ratio A (defined below). 

|| The roughness exponents \ u and \b of the fields u and b are related to Xh and X<t> : Xh = Xu + lj X<t> = 
Xb + 1- 
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4-1. Review of the KPZ Equation 

The KPZ equation for surface growth is one of the simplest non-linear generalization 
of the diffusion equation and serves as a paradigm for phase transitions and scaling in 
non-equilibrium systems; see for example Refs. [UJ,[26] for extensive reviews. Our GBM, 
Eq. (JHJ), reduces to the KPZ equation for = 



The field fr(x, t) physically represents the height profile of a growing surface. The 
KPZ equation has by now been studied by a broad variety of approaches. These include 
dynamic RG [27], Monte-Carlo simulations of the equivalent lattice-gas models [21], and 
by mapping onto the equilibrium problem of a directed polymer in a random medium 
[T7j . The main results concerning the statistical properties include 

• Due to the Galilean invariance the scaling exponents xh and z characterizing the 
rough phase follow an exact relation Xh + z = 2 [T5], [23] . 



• The scaling exponents for the rough phase in Id are exactly given by Xh — 1/2, z — 
3/2 [UJ [231 [26] . 



• For spatial dimensions d = 2 + e with e > there is a phase transition from 
a smooth to a rough phase. The exponents in the smooth phase are exactly 
given by Xh = (2 — d)/2 and z = 2 [28]. The values of the exponents in 
the rough phase are controversial. Functional renormalization group studies [T7J 
and equivalent mode coupling analyses in terms of a small- x expansion [23] give 
Xh — (4 — d)/Q and z = (8 + d)/Q. However, a recent critical analysis of 
the mode coupling equations by Canet and Moore [29] modifies these values. 
They obtain z = 2 - (4 - d)/A + 0((4 - cf) 2 ), Xh = 2 - z near d = 4 and 
z = 4/3 + d/3 + 0(d 2 ), Xh = 2 — z near d — 0. These scaling exponents describe 
the rough phase at d = 1 and 2 < d < 4. They further obtain a new set of scaling 
solutions of their mode-coupling equations for d < 2. 

The GBM discussed in this work are expected to exhibit much richer behavior, since 
in addition to advection and diffusion it contains a feedback term (the term A2V6 2 ) 
in Eq. (Q. Furthermore, although equations (Q and ffTUj) are invariant under parity 
inversion, since the fields u and b have different properties under parity inversion, an 
intriguing possibility of breaking parity in the statistical steady state by the presence 
of a nonzero cross-correlations of u and b (created by suitable choices external forces) 
exists. We discuss some of these issues below. 

4-2. Dynamic renormalization group studies in d-dimension 

In this section we employ dynamic renormalization group (DRG) methods to understand 
the long-time and large-scale physics of the model equations ([9]) and (flOi) . Before going 
into the details of our calculations and results we elucidate the continuous symmetries 
under which the equations of motion remain invariant. As shown in appendix 19.11 these 



dh 
~dt 



2 h. 



(32) 



Scaling and universality in coupled driven diffusive models 



12 



allow us to construct exact relations between different vertex functions which in turn 
impose strict conditions on the renormalization of different parameters in the model. 
Here we list the symmetries and summarize the consequences on the renormalization of 
the model. 

• The model shows Galilean invariance when Ai = A3 = A, i.e., the equations of 
motion are invariant under the continuous transformations 

u'(x, t) = u(x + Xu t, t) , 

b'(x, t) = b(x + Auot, t) . (33) 

This invariance implies that the coupling constants Ai and A3 do not renormalize 
in the long wavelength limit. 

• There is a rescaling invariance of the field 

b->Ab, A 2 ^ A 2 /VI. (34) 

This ensures that also the coupling constant A2 does not renormalize in the long 
wavelength limit. This can be formulated in a more formal language (see Appendix 

EH). 

Summarizing, none of the coupling constants Ai, A2, A3 renormalizes. Hence, 
in perturbative renormalization group treatments, if we were to carry out a 
renormalization-group transformation by integrating out a shell of modes Ae~ l < q < A, 
and reseating t — > e lz t, bi — > e xl bi, the couplings Ai, A 2 and A 3 

would be affected only by naive rescaling. Thus, rescaling of space and time can be 
done in such a way as to keep the coupling strengths Ai, A 2 , A3 constant. 

We perform a one-loop dynamic renormalization group (DRG) transformation on 
the model Eqs. ([9]) and ( fTUl) with the correlations of the noise sources specified in 
Eqs. (1171) . (1181) and ( [191 ). The cross-correlation function is imaginary and odd in 
wavevector k: It is proportional to sgn(k) and hence is non-analytic at k = 0. Since 
perturbative expansions as in DRG analyses used here are always analytic at k = 0, non- 
analytic terms of the form sgn(h) are not generated. Hence there are no perturbative 
corrections to the cross-correlations. Furthermore, the model equations ([9]) and (TTOl) 
and the variances of the noise sources (fTTl) , (ITSl) are invariant under inversion (reversal 
of parity), rotation and the interchange of Cartesian co-ordinates i and j. The last two 
invariances, under rotation and interchange between % and j respectively, are broken only 
by the choice of the noise cross-correlations (fl9~i) . of the external stochastic forces which 
can be controlled from outside separately. This is reflected in the fact the presence of 
the symmetric (anti-symmetric) noise cross-correlations does not lead to the generation 
of the anti-symmetric (symmetric) noise cross-correlations. This allows us to explore the 
effects of symmetric and anti-symmetric noise cross-correlations separately. It should 
be noted that the calculations presented here are done at a fixed dimension d, instead 
of as an expansion about any critical dimension. 
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Symmetric cross- correlations. We first consider the case when the (bare) noise cross- 
correlations are fully symmetric (finite D®), i.e., no anti-symmetric cross-correlations 
are present (D® = 0). We perform a renormalization group transformation as outlined 
above. The resulting RG flow equations are presented in terms of the renormalized and 
rescaled variables 

(u, n) - (u, fi)e^ 1 , (D u , D b , D s ) - (D u , D b , D s )e^ d - 2 ^ 1 . (35) 

We obtain the following differential flow equations for the running parameters 
1/(0, MO, A»(0, D bV) (with Ai = A 2 = A 3 = A) 



du 

~dl 

dji 
~dl 
dD u 
dl 

dD b 
~df 
d\ 
111 



dG 



l' 



D, 



D h 



z-2 



d 
G2 



4 

- d 



■2 Xh 
2 Xh - 



d 
d- 



(1 + AI*) 

3 4~ Pm 
1 + P m ? 

1 + A 2 P n 



G 
7 



A(l + 3P„ 

(i+p m y 

4N«P 2 



p 

J. rr 



d + G 



P 2 



4 p3 



p 



P m f A 



2]- 



(36) 
(37) 
(38) 
(39) 
(40) 



The parameter D s does not receive any fluctuation corrections and is affected only by 



naive rescaling. Here, we have introduced an effective coupling constant G 



X 2 Du 



and 



two amplitude ratios, A = ^ and N s = (-^p) 2 , characterizing the relative magnitude of 
the noise amplitudes for the magnetic field and the (symmetric) cross-correlations with 
respect to the noise amplitude of the velocity field, respectively. Further, P m = u/fj, is 
the renormalized Prandtl number. We find, from the flow equations ( 1361) and ( 1371) . at 
the RG fixed point v — fi, i.e., we have P m = 1 at the RG fixed point, regardless of the 
values of the bare viscosities. Henceforth, we put v = fi in our calculations below. Flow 
equations for the effective coupling constant G and the amplitude-ratio A = D b /D u may 
be obtained from its definition above and by using the flow equations ( 13611401) . They are 



dA 

~dT 

dG 

~df- G 



A[(1 + A 2 + 2N S )-2(1-N S /A)] 



G 



2-d + 2G 



2d -3 
2d ' 



iV, 



(41) 
(42) 



At the RG fixed point renormalized parameters are scale invariant (i.e., do not 
receive fluctuation corrections anymore under further mode eliminations); we then set 
the LHS of Eqs. (I36H40I) to zero. These yield (a * denotes fixed point values), 

2de , 



CT e < 



2d- 3 v 

to the lowest order in N s with e 
G* = (1 + N s ) we find 

A* = 1 - 2iV e . 



(43) 



d — 2. When G* = 0, A is undetermined and when 



(44) 
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Note that the fixed point value of the effective coupling constants G and the amplitude 
ratio A explicitly depend on the strength of the noise cross-correlations N s . We show 
below that the parameter N s is marginal at the RG fixed point, i.e., N s can have variable 
values at the fixed point. In the rough phase at d — 1, we find G* = 2(1 + N s ) 
as the stable fixed point and an amplitude-ratio A* = Db/D u = 1 — 2N S . This 
implies that the non-linearities are relevant and the asymptotic scaling properties of 
the correlation functions are different from those for the corresponding linear model. 
For d = 2 + e with e > 0, we obtain G* = as a stable and G* = 4e(l + N s ) as 
an unstable fixed point indicating a smooth-to-rough transition. In the smooth phase 
(G* = 0) the nonlinearities are irrelevant, the scaling properties are determined by the 
corresponding linear equations and we find dA/dl = 0, i.e., A does not change under 
mode elimination and is simply given by D®/D®, the bare amplitude ratio. At the 
roughening transition, G* = 4e(l + A,), the amplitude ratio becomes A* = 1 — 2N S . 
Also note that the value of the coupling constant G* at the critical point for d = 2 + e 
increases with increasing A,, i.e., with increasing symmetric noise cross-correlations. 
These, therefore, suggest that the presence of symmetric noise cross-correlations helps 
to stabilize the smooth surface against roughening perturbations. For d > 2 beyond the 
critical point (roughening transition point) there is presumably a rough phase which is 
not accessible by perturbative RG. This is reminiscent of the analogous problem in the 
KPZ Equation [23J. Further we obtain \h — 1/2, z = 3/2 in the rough phase at d — 1 
and x = -C(e 2 ), z = 2 + 0(e) 2 at the roughening transition for d = 2 + e. Scaling 
exponents for the rough phase cannot be obtained by perturbative RG. 



Anti- symmetric cross-correlations . Having discussed the effects of symmetric cross- 
correlations, we now proceed to analyze the effects of the anti-symmetric cross- 
correlations in a DRG framework. Therefore, we now have a finite D® in the bare 
noise cross-correlations with D® being set to zero. In this situation, fields u and b, and 
forces f and g are no longer expressible as gradients of scalars. Hence, Equations (J9j) 
and ( fTOl) cannot be reduced to ( fl~4"|) and ( |T5l) . Therefore, we work with Equations (jHJ) 
and (TlOl) directly. We follow the same scheme of calculations as above. The resulting 
RG flow equations are 



dv 
~dl 
dfi 
~dl 
dX 

111 
dD 
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(45) 
(46) 
(47) 
(48) 

(49) 
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The parameter D aj like D s above, does not receive any fluctuation correction and is 
affected only by naive rescaling. Note that the flow Eqs. (145 p and (I46p have the same 
form as Eqs. ( |36|) and ( 1371) . the corresponding flow equations in the symmetric cross- 
correlations case, except that we now denote the new dimensionless coupling constant by 
G a = \ 2 D u /is 3 . Here, N a = (D a /D u ) 2 . Further, as in the symmetric cross-correlations 
case, P m = 1 at the RG fixed point. However, the flow equations (jUJ) and fj4~9]) are 
not identical to their counterparts (1381) and (|39|) in our discussions on symmetric cross- 
correlations above: the fluctuation corrections contributing to D u and Db arising from 
the anti-symmetric cross-correlations have the same sign in this case unlike the case 
with symmetric cross-correlations. We find that the parameter A = Dt>/D u is unity 
for all values of D a , i.e., A = 1 at the RG fixed point. Further, one can obtain a flow 
equation for the the coupling constant G a = X 2 D u /u 3 with the help of Eqs. (1451 - 1471) . It 
is 



dG a _ 

~dT- Ga 



2d — 3 

2-d + 2G„| — — + N a 



2d ■-■ ■ < 5 °» 

At the RG fixed point dGJdl = yielding G a = or G a = e + 0(N a )] . 
The value G a = corresponds to the smooth phase, as in the KPZ case, whereas 
G* a = e [gjzjj + 0(N a )\ is an unstable fixed point indicating a smooth-to-rough phase 
transition. In the smooth phase the scaling exponents z = 2, \ = (2 — d)/2 which are 
identical to their KPZ counterparts. Further at the phase transition point the exponents 
are independent of N a : z = 2 + 0(e) 2 , x = 0(e) 2 ; again the statistical properties of the 
rough phase cannot be explored by perturbative RG. 

Note that our above conclusions on obtaining continuously varying amplitude-ratios 
in the rough phase in Id and at the smooth-to-rough transition for d > 2 for finite 
symmetric cross-correlations rest on the requirement that N s and hence D s can have 
variable values at the RG fixed point. This can happen if D s is marginal at the RG 
fixed point. Therefore, to complete our analysis we now proceed to demonstrate that 
the parameter N s is strictly marginal at the RG fixed point, even beyond linearized RG. 
In our notations D s (k.)/u is the amplitude of C-j(k,t), the symmetric part of the cross- 
correlation function matrix, with D s (k)D s (k) = 4D 2 . Since C^(k,t) is an odd function 
of the wavevector k, we must have, for consistency, C^-(k = 0,t) = 0. Therefore, there 
must be a length scale l c (which itself diverges in the thermodynamic limit) such that 

C« (k, t = 0) = 2iD s i] (k)k~ d ~ 2x /v for k -> (51) 

up to a scale I, 1/k ~ I < l c — > oo, and 

C* (k, t = 0) = 0, at k = (l c < I = oo). (52) 

Under rescaling we have D u (l) = D u l z ~ 2x ~ 3 as / — > oo. In contrast, under the same 
rescaling cross- correlation C? (k, 0) = 2iD s {k)k~ d ~ 2x kikj/v if k' = k/b > l~ l but 
is zero if k! < Thus the true scaling regime is / < l c — >• oo and at I < l c , 

D s (l) ~ D s l z c - 2x ~ d - 2 . The latter does not receive any fluctuation corrections under 
mode elimination, as we argued above, and hence is arbitrary because it depends on 
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D s = D S {1 = 0) and therefore marginal. Therefore, N s = (^) 2 is also marginal. The 
marginality of N a = (^ £L ) 2 can be argued similarly. We close this Section by summarizing 
our results obtained from the DRG scheme: 

• We obtain the amplitude-ratio D b /D u = 1 — 2N S in the presence of finite symmetric 
cross-correlations N s in the rough phase at d = 1 and smooth-to-rough transition 
at d — 2 + e. The values of the scaling exponents are unchanged from their values 
for the KPZ equation. 

• In the presence of the anti-symmetric cross-correlations, the amplitude-ratio 
D b /D u = 1 at d = 2 + e. The scaling exponents are unchanged from their KPZ 
values. 

Clearly, the DRG scheme, as in the KPZ equation, fails to yield any result concerning 
the strong coupling rough phase at d = 2 + e. We now resort to the SCMC and the 
FRG schemes to obtain results for the rough phases at d = 2 + e. 

4-3. Self- consistent mode- coupling analysis in d-dimensions 

In a self-consistent mode-coupling (SCMC) scheme perturbation theories are formulated 
in terms of the response and correlation functions of the fields u(x, t) and b(x, t). They 
are conveniently expressed in terms of self-energies and generalized kinetic coefficients. 
As before, without any loss of generality we assume v = fi, i.e., the magnetic Prandtl 
number P m = 1. This guarantees that there is only one response function and it can be 
written as 

G-l(k,u) = -iu-Z(k,uj), (53) 
The correlation functions are of the form 

C% b (k,u) = 2D Uib (k,u)k l k J \G(k,u)\ 2 , 

Cg(k jW ) =2iD x (k,u)ij\G{k,u)\ 2 , (54) 

where (k, uj) stands for the cross correlation functions of u and b. In the scaling limit, 
in terms of wavevector k and frequency u the self-energy and the correlation functions 
exhibit scaling forms characterized by the scaling exponents z and x an d appropriate 
scaling functions: 

Cf{k,u) = 2D u , b k i k j k- 2 x- 2 - d - z f u , b (u;/k z ), 

Cg(k,w) = 2zD x (k)k i k J k-^- 2 - d - z f x (u/k z ). (55) 

In diagrammatic language a lowest order mode-coupling theory is equivalent to a self- 
consistent one-loop theory. The ensuing coupled set of integral equations are compatible 
with the scaling forms above. To solve this set of coupled integral equations we 
follow Ref. [23] and employ a small-x expansion. This essentially requires matching of 
correlation functions and the self-energy at zero frequency with their respective one-loop 
expressions. We consider the following two cases separately: 
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Symmetric cross- correlations:- We consider the case when C^(k, u;) = Cj*(k, u). This 
implies that the fields u and b can be expressed as gradients of scalars: u = V/i and 
b = V0; note that is actually a pseudo scalar. Such choices imply that the fields u 
and b are irrotational vectors. The corresponding equations of motion in terms of the 
fields h and (f> are given by ((Hj) and ( fl"5l) . 

In this case the zero frequency expressions for the correlators and the response 
function become 



£(k,0) 



Tk z 



C£ 6 (k,0) = 2 



D 



u.b 

ij h h h 



-d-2 X -z 



Q(k,0) 
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D s (k) 



-d-2x- 



(56) 



where C?-(k, 0) is the symmetric part of the cross correlation function of u and b. We 
also define D s (k)D s (k) = (D s ) . In an SCMC approach vertex corrections are neglected 
which are exact statements for the present problem in the zero-wavevector limit. Lack 
of vertex renormalizations in the zero-wavevector limit yield the exact relation between 
the scaling exponents x an d z, as in the case of the noisy Burgers/Kardar-Parisi- Zhang 
equation [23]. In the context of the Burgers equation in 1+1 dimension Frey et al 
showed |25] , by using nonrenormalization of the advective nonlinearity and second order 
perturbation theories that the effects of the vertex corrections at finite wavevectors on 
the correlation functions are small. Presumably the same conclusion regarding the effects 
of vertex renormalization at finite wavevectors follows for this model in the present 
problem also. However, a rigorous calculation is still lacking. With these definitions the 
one-loop self-consistent equations yield the following relations between the amplitudes. 
For the self-energy we obtain (without any loss of generality we set Ai = A2 = A3 = A) 

a^ = (27p2rf (1 + irj' (57) 

and for the one-loop correlation functions 

r 2 1 s d 1 
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Here Sd is the surface of a d-dimensional sphere. From Eqs. (l57|) and (|58|) we obtain 



(5f 



1 = 0. 



(59) 



where N s = (D s /D u ) 2 is a dimensionless ratio as defined above. Since the ratio D b /D u is 
positive semi-definite, in Eq.( l59l the range of N s is determined by the range of positive 
values for D b jD u starting from unity (obtained when N s = 0). Thus for small N s we 
can expand around zero and look for solutions of the form A = D b /D u = 1 + aN s , 
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such that for N s = we recover D u = Dj, (the result of Ref-jTH]). We obtain a = —2, 
i.e., Db/D u = 1 — 2N S , implying that within this leading order calculation N s cannot 
exceed 1/2, i.e., D s < D u /y/2. An important consequence of this calculation is that 
the amplitude ratio Di,/D u is no longer fixed to unity but can vary continuously with 
the strength of the noise cross-correlation (renormalized) amplitude D s . Our results 
from this section are in agreement and complementary to the those obtained in a DRG 
framework above (see Sec. 14.21) . These results are already confirmed by a one- loop DRG 
calculation for the rough phase at d — 1 (see Sec J4. 2ft . In addition, the application of one- 
loop DRG demonstrates that the above results are valid at the roughening transitions 
to lowest order in N s in a d = 2 + e expansion as well. 

In contrast, the scaling exponents x an d z are not affected by the presence of cross 
correlations. From the above one-loop mode-coupling Eqs. (l57|) and ( |58l) we obtain 
X = — 1/2 and z — 3/2 in d — 1 dimensions from our SCMC which are same as obtained 
by DRG calculations. Further equations (1371) and (1581) yield the following values for the 
scaling exponents in the strong coupling regime which are unaffected by the presence of 
symmetric cross-correlations. 

1 d 4 d , \ 

X = ; z = - + -. 60) 

A 3 6' 3 6 V ; 

These are identical to those obtained by Bhattacharjee in a small- x expansion [23] and 
it is still controversial whether these values for the exponents actually correspond to 
the usual strong coupling case. Recently, Canet et al performed a more critical analysis 
of the self-consistent mode-coupling equations for the KPZ Equation [2H] and showed 
the corresponding mode-coupling equations have two branches (or universality classes) 
of the solutions: the F branch having the upper critical dimension d c = 4 and the S 
solution with d c — 2. The F solution is believed to correspond to the usual rough phase 
and the S solution has been discussed in some calculations on the directed polymer 
problem [30]. Our solutions or rather Bhattacharjee's small x expansion yields d c = 4 
as the F solution and agrees with the F solution at d = and d = 1 as well. At other 
dimensions there are small quantitative differences between the values for z. 

Antisymmetric cross- correlations: So far we have restricted ourselves to the case where 
the vector fields u and b are irrotational. If however the fields a = u, b are rotational 
and have the form 

a=VxV a + VS a , (61) 

with vectors V a being cross-correlated but the scalars S a uncorrelated then the variance 
D*j(k) satisfies 

^(k) = -£>*(-k) = ^(k) = -[L>*.(k)]*. (62) 

This is the antisymmetric part of the cross-correlations. Choices (1611) ensures that the 
vectors u and b are no longer irrotational. The corresponding (renormalized) noise 
strength D a is formally defined through the relation 

D*Jk)D*(-k)=4Dlk\ (63) 
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Similar to the previous case, in the scaling limit (zero frequency limit) the self energy 
reads E(A;, u = 0) = Tk z , the correlation functions are C^(k,uj = 0) = kikjD u k~ d ~ 2x ~ z , 
C^(k,uj = 0) = kikjDbk~ d ~ 2x ~ z , and the antisymmetric part of the cross-correlation 
function reads C?Ak,u = 0) = D^(k)k 



-ax- 



Following the method outlined above we obtain 
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(65) 



Equations (|64"|) and fl65|) give D u /Dj, = 1 at the fixed point for arbitrary values of 
N a = (D a /D u ) 2 . Hence no restrictions on N a arises from that. In contrast to the effects 
of the symmetric cross-correlations, the exponents now depend continuously on N a . To 
obtain the scaling exponents we use that D b /D u = 1 and equate Eqs. (l6"4l) and (16"5"|) . To 
leading order, we get 

1 d 



N a d 



N a d 



y = 1 , z = — I . (66) 

A 3 6 6 ' 3 66 K J 

These exponents presumably describe the rough phase above d > 2, with the same 

caveats as above [28]. With increasing D a the exponent \ grows (and z decreases). 

Obviously this cannot happen indefinitely. We estimate the upper limit of N a in 

the following way: Notice that the Eqs.Q and fflUj) along with the prescribed noise 

correlations (i.e., equivalently the dynamic generating functional) are of conservation 

law form, i.e. they vanish as k — ► 0. Thus there is no information of any infrared cut off 

in the dynamic generating functional. Moreover, we know the solutions of the equations 

exactly if we drop the non-linear terms (and hence, the exponents: x = l — d/2, z = 2). 

Therefore, physically relevant quantities like the total energies of the fields u and b fields 

^1, f, (u(k, £)u(— k, £)) and J fc (b(k, t)b(— k, £)), remain finite as the system size diverges, 

and are thus independent of the system size: In particular 



(u(k, *)u(-k,t)) 

;b(k,t)b(-k,t)) 



d d k k 2 ~ d+2x , 
d d k k 2 ~ d+2x 



(67) 



which, for x — —d/2 (the exact value of x without the nonlinear terms), are finite in 
the infinite system size limit. Since the non-linear terms are of the conservation law 
form, inclusion of them cannot bring a system size dependence on the values of the 

These are kinetic and magnetic energies when u and b are interpreted as Burgers velocity and 
Burgers magnetic fields. 



Scaling and universality in coupled driven diffusive models 



20 



total energies. However, if x continues to increase with D a at some stage these energies 
would start to depend on the system size which is unphysical [3T] : 



Therefore, in order to make our model meaningful in the presence of anti-symmetric 
cross-correlations, we have to restrict D a to values smaller than the maximum value 
for which these energy integrals are just system-size independent: This gives iV™ ax = 
|(c//2 + 1). Note that the limits on N s and N a impose consistency conditions on the 
ratios of the amplitudes of the measured correlation functions but not on the bare noise 
correlators. We can use the values of the dynamic exponent to estimate the upper critical 
dimension d c of the model in the presence of the antisymmetric cross-correlations. From 
our expressions ( l66l) the dynamic exponent z increases with the spatial dimension d for 
a given strength of the antisymmetric cross-correlation N a . The value of d c is given 
by the dimension in which the dynamic exponent z attains a value 2 equal to its value 
without the nonlinear term. Clearly, from expressions (1661) z = 2 yields d c = 4/(1 — N a ). 
Therefore, the antisymmetric cross-correlations have the effects of increasing the upper 
critical dimension of the model. 

Antisymmetric cross-correlations stabilize the short-range fixed point with respect 
to perturbations by long-range noise with correlations of the form (in the Fourier space) 
k~ y , y > 0. This can easily be seen: In presence of noise correlations sufficiently singular 
in the infra-red limit, i.e. large enough y, the dynamic exponent is known exactly [28| [32] : 
For a sufficiently large y the one-loop corrections to the correlators scale same as the 
bare correlators for zero external frequency; the one-loop diagrams are finite for finite 
external frequencies. Thus they are neglected and this, together with the Ward identities 
discussed above yield z\ x = — |. The short range fixed point remains stable as long 
as z ST < z\ T which gives y < —2 + (1 + N a )d/2. Hence we conclude that in the presence of 
antisymmetric cross-correlations a long range noise must be more singular for the short 
range noise fixed point to loose its stability or in other words, antisymmetric cross- 
correlations increases the stability of the short range noise fixed point with respect to 
perturbations from long range noise sources. 

We close this section by summarizing our results obtained from the SCMC 
calculations: 

• The SCMC method yields results about the rough phase at d = 1 and d > 2. 

• We find that the amplitude Db/D u decreases monotonically from unity as the 
symmetric cross-correlations, parametrized by N s increases from zero. Our result 
here is in agreement with that obtained from the DRG method for d = 1 . The ratio 
Db/D u is unaffected by the antisymmetric cross-correlations. 

• Our SCMC calculations yield for the scaling exponent also. In the presence of the 
antisymmetric cross-correlations parameterized by N a they are: x = ~ f — f + 




(68) 
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PP , z = | + | — pp. The scaling exponents are unaffected by the symmetric 
cross-correlations. 

• The maximum value of the parameter N s is obtained by setting Db/D u to zero, 
while the maximum value of N a is obtained by setting \ t° zero i n an y dimension. 
The minimum values for both of them are zero. 



4-4- Functional renormalization group analyses on the model 

In this Section we study the model Eqs. (Q and ([TO]) in a functional renormalization 
group (FRG) framework. This study is complementary to our DRG and SCMC studies 
above. In Section [3] it has been discussed that the model Eqs. © and (fTPl) . for the 
bare Prandtl number = i^/Vo = 1, reduces to two KPZ equations [see, Eqs. (12~T]) for 
their Id representations] representing two growing surfaces hi(x.,t) and h 2 (x, t) which 
are coupled by noise sources [see, Eqs. fl26l) for the noise sources in Id]. Such equations 
in general ^-dimensions are 



d t h + l/2(V/n) 2 = z/ V 2 /ii + ^i, 

d t h 2 + l/2(V/i 2 ) 2 = z/ V 2 /i 2 + V>2- (69) 
Here, the noise correlations in arbitrary dimension d are given by 

<^i(k,^i(-k,0))=2A)5(0, 
<V 2 (M)V2(-k,0)> =2D 5(t), 

(^i(k, t)^ 2 (-k, 0)) = 2A,5(t) + 2iA)(k)£(i), (70) 
omitting a formally divergent factor = 0) in each of the variances above. Note that 
in the third equation of (ITDjl the cross-correlation of the noise sources ipi and ^ 2 has a 
real and an imaginary parts, where as in its one- dimensional version, used to introduce 
our lattice-gas models, given by Eqs. (I24H26I) the cross-correlation has no real part. This 
is because in the bare theory even if the real part of the noise cross-correlation is zero 
it would be rendered non-zero self-consistently in the presence of the imaginary part. 
In other words, the imaginary part gives rise to the real part in the (one-loop) self- 
consistent theory. The real part, however, remains zero self-consistently if there in no 
imaginary part in the bare noise cross-correlations. 

We begin by applying the well-known Cole-Hopf transformation [33] to the Eqs. 
(I69p : /ii,2(x, t) = (2v/X) hi Zip. These transformations reduce the Eqs. ( |69l) to 

dtZ x = v^ 2 Z x + 

d t Z 2 = u V 2 Z 2 + V 2 Z 2 , (71) 

where V\ = V 2 = j^. Equations (ITT]) can be interpreted as the equations for the 
partition functions Z\ and Z 2 for two identical directed polymers (DP), each having d 
transverse components, in a random medium whose combined Hamiltonian is given by 



H= dt 



2 V dt J 2 [ dt 



V r i(xi,t) + V r 2 (x 2 ,t) 



(72) 
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Functions Vi(x 1; t) and V2(x 2 , t), for this generalized DP problem, are to be interpreted 
as quenched random potentials experienced by the two DPs embedded in them. Clearly, 
the potentials V\ and V2 are Gaussian distributed with zero-mean and variances given 
by Eqs. (1701) . The coordinate t in the Hamiltonian H in Eq. (1721) . which denotes the 
physical time in the coupled surface growth problem, now becomes the arc-length of 
the DPs; xi(i) and x 2 (t) are the transverse spatial coordinates of the two DPs. In the 
Hamiltonian (I72|) the first two terms are the energies of the two DPs due to transverse 
fluctuations (elastic energies) which are minimized if the DPs are straight, and Vi, V2 
are the potential energies due to the quenched disorder which can be minimized if 
the DPs follow the minima of the potential landscapes (and hence they will not be 
straight). Thus there will be competition between the two opposite tendencies and 
there will be different phases depending upon which one wins in the thermodynamic 
limit. Due to the structure of the noise correlations given by the Eqs. ( 1701) it is clear 
that the cross-correlations of the quenched random potentials V\ and V 2 have a part odd 
in wavevector k, suggesting that the disorder distribution of the embedding disordered 
medium lacks reflection symmetry, i.e., it has a chiral nature. The phase diagram of 
two DPs in a reflection-symmetric random environment has been discussed in Ref.[18j. 
In the present work we include the effects of chirality and discuss its consequence on 
the statistical properties of the two DPs. Thus, our studies of the Eqs. (Q and (fTUl) 
can equivalently, in terms of the Directed Polymer (DP) language, be considered as 
investigating the phase diagram of the following toy model: Let us assume that the two 
DPs A and B are embedded in a random medium which has two kinds of pins A and 
B which pin polymers A and B respectively. Both the pins are distributed randomly 
with specified distributions. Furthermore, the pins A and B may have some correlations 
in their distributions, or may not have. If they do have, then, the effects of such 
correlations should be modeled by the cross-correlations of the type we are discussing 
here. Physically, if there is a pin A somewhere, then a positive correlation between the 
distribution of pins A and B would indicate that a pin B is likely to be found nearby. 
Since pins are the places where polymers are likely to get stuck, then according to the 
above, if polymer A is stuck somewhere, then polymer B is also likely to get stuck nearby 
with a probability which is higher than if pins A and B had no correlations among their 
distributions. This, in some sense, creates an effective interaction between polymers A 
and B (since they are more likely to get pinned at nearby places). We elucidate the 
resulting effects in an FRG framework. 

In the present problem, since the two DPs are not interacting with each other 
directly, the total partition function Z of the two DPs, for a given realization of the 
pinning potentials, is then given by the product of the individual partition functions: 



Here, T = l/u is the temperature. Therefore, following the standard replica method 
[M] the free energy of the system, after averaging over the distribution of the random 




(73) 
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potentials V\ and V 2 , is given by 

(F) = (InZ) = lim N ^ - 



(Z N ) 



N 



(74) 



In order to facilitate the usage of the standard FRG method we generalize the 
correlations of the random potentials (ITUI) . similar to the corresponding FRG treatment 
for the single DP problem [T7j, in the following way: 

(V 1 {x,t)V 1 (Q,Q))=R 1 (x)S(t), 

(y 2 (x,t)y 2 (o,o)) = i2 2 (x) ( j(t), 

(Vi(x, t)V 2 (0, 0)> = R(x)5(t) + R x (x)5(t). (75) 

Here, following Ref. [17], the spatial (^-functions in (1701) have been replaced by short 
range functions R±, R 2 , R for convenience. The function i? x (x) is an odd function of 
x representing the imaginary part of the cross-correlations (^i(k, t)%p 2 {— k, 0)). Hence, 
-Rx(O) = 0. Note that in the DP-language functions -Ri(x), -R 2 (x), -R(x) and -R x ( x ) 
are proportional to the noise variances in Eqs. ( ITUI) and hence to the corresponding 
correlators. With the above definitions and notations, we have 



IL„ „Vx lfJ Vx 2 v exp 



1 

2T 



dtY, r 



c?x 2q . 



x exp 



x exp 



— S Q)( gi2i(x lc 



X1/3) + — S Qi/ g^ 2 (x2/3 



1 1 

^2-S Q]j a_R(x lQ — x 2/ g) + — H a pR x (x lc 



" X 2/3) 



(76) 



Here, indices correspond to the replica indices arising out of the replica method 
used above, representing identical copies of the same system. Clearly, averaging over 
the distribution of the potentials lead to generation of terms with mixed replica indices 
- systems having different replica indices now interact with each other. 

In order to set up the functional renormalization group (FRG) calculation for 
establishing the long wavelength forms of the disorder correlators we rescale x x 2 — ► 



e'xi 2 and t e & such that t 



l x i, 2 |^, relating longitudinal and transverse 



fluctuations. Clearly, the exponent £ is the inverse of the dynamic exponent z: ( = 1/z. 
Such a rescaling yields for temperature T — > e( 1_2< ^T. The differential flow equation 
for T then reads 

^ = (1 - 2C)T. (77) 

Hence, if the disorder induced roughening dominates over thermal roughening (( > 1/2) 
we have T — >• under renormalization and the long wavelength physics is governed by a 
zero-temperature fixed point [TT| . In a functional renormalization group (FRG) analysis 
one splits the degrees of freedom (here xi and x 2 ) into their long and short wavelength 
parts: 

x 2 |+x 2 >. (78) 



Xl„ + x i^, x 2/3 
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We write the degrees of freedoms x x and x 2 as in (1781) and consider the part of H 
quadratic in the short wavelength parts of the degrees of freedom {H > ). We expand the 
disorder potential terms containing x^ , x 2 in the exponential of ( 1761) up to the second 
order [IT], average over the x^ , x 2 , and neglect terms containing three- replica indices 
due to their irrelevance p2] to obtain 

H> = T 1 + T 2 + T 3 , 



(79) 



where 



4T 2 J g 4^ a >V 



R'((^< - Xl <) 2 - 2S Q , AtJ R' 1 / (0) J R / 1 '(xi< - x x 
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2^J 



(80) 

with g being the Fourier conjugate variable of t. Note that arguments of Ti, T 2 and 
T 3 are, pure xiq, — x^, x 2q , — x 2ai and mixed xi a — x 2m respectively. Comparing with 
the existing (bare) terms in H < we find that Ti, T 2 and T 3 , renormalize, respectively 
-Rf(xiQ, — x lM ), i? 2 < (x 2Q , — x 2m ) and _R < (x lQ , — x 2/ J. Since the all of Ti, T 2 and T 3 are 
even under inversion of their arguments, there are no corrections to R x . Corrections 
Ti and T 2 , to Ri and i? 2 respectively, have identical forms and are same with the 
corresponding corrections in the single DP case [IT]. This is expected, since before 
disorder averaging, the free energies of each of the DPs, like the free energy of the 
single DP problem, follow the usual KPZ equation. In the expansion of the disorder 
correlation terms in the exponential of ( 1761) the terms in the first order of the expansion 
do not contribute. This is because the above expansion is essentially perturbative in 
T as in the single DP case [TTj : T flows to zero under renormalization [see eq. (1771)]. 
Hence the first order terms having an uncompensated power of T flow to zero and are 
irrelevant (in an RG sense). This feature is same as in the single DP case [TTJ. Note 
that we have made use of the fact that R x {0) = while arriving at the expression (179]) . 
Different terms in (1791) contributes to the fluctuation corrections to Rf, Rf, R < which 
can be identified by their arguments. Note that there are no corrections to R x which 
is reminiscent of the lack of fluctuation corrections to the noise cross-correlations in the 
Eqs. ([9]) and ( flOl) in the long wavelength limit. 
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In the next step, we argue that all of the functions R±, R 2 , R, R x are characterized 
by the same scaling behavior in the long wavelength limit. This is because all of them 
are proportional to various noise variances in the model given by Equations ( TT7T) . ( {TBI) 
and (fT9l) . Now all the correlation functions in the model Eqs. (Q and (flOl) . in stochastic 
Langevin descriptions, are proportional to the noise variances ( IT71) . (fTBl) and ( {T9lh 
Further, these correlation functions, due to the symmetries of the GBM model, have the 
same scaling behavior in the hydrodynamic limit, characterized by a single roughness 
(x) and dynamic (z) exponents. Hence, the effective noise variances, and therefore, the 
functions Ri, R 2 , R and R must have the same scaling behavior in the long wavelength 
limit. With the rescaling of t and x mentioned above and the i?-functions scale as 

R -»■ [1 + (3 - ±()5l]R, (81) 

where R in the above stands for all of R±, R 2 , R, R. These then yield the following 
differential flow equations: 



dR x 



dR 2 
dl 



dR 



dR, 
dl 



(3 - 4C)i?i + (xR[ 
R'{{x) 



1 



(R'lY - R"(0)R'l 



n — 1 



R[ 



x 



in 



i)R';(oy 



X 



(3 - 4()R 2 + (xR' 2 

1)R'M- 
4()R + (xR' 



(R' 2 y-R' 2 \0)R' 2 > + 



n 



Rn 



r 



n 2 



[n 



X 



(3 



+ 



R' 



x 



— [n 



-{R"f 



l)i?"(0) 



R"(x) 



n 



x 



R 

x 



{3-A()R x + (xR' x . 



(82) 



In the Eqs. (152^) above ""' denotes a derivative with respect to x, the argument of the 
functions R± etc. Note that the functional flow equations for the functions R\ and R 2 
are identical to each other which is expected on the ground of symmetry between the 
equations of motion fl69|) or (ITT|) . At the RG fixed point all the partial derivatives with 
respect to the scale factor / is zero yielding 



(3-4C)i?i 



(xR' 1 + ~(R'tf-R';(0)R>; 
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- (n - 1)^(0) 
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In order to proceed further, we make the following choice without any loss of generality: 
R(x) = ^R\(x) = , ~fR2(x) in the long wavelength limit, where 7 is a numerical constant. 
We further choose R x (x) 2 = T Ri(x) 2 = T R 2 (x) 2 in the long wavelength limit, where 
To ia another numerical constant. These parametrizations are consequences of the 
symmetries of the GBM model, which ensure, as we have argued above, functions 
Ri, R2, R and R are proportional in the thermodynamic limit. We determine below 
a relation between 7 and r .In terms of the parameters 7 and T , then, the first and 
the third in the Eqs. (1821) at the fixed point reduce to 



(3 - 4C)i2i 
(3 - K)Ri 



Since the equations in 



7 + r = 7 



Ri + (xR[ + 
Ri + (xR[ + 
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R'i = 0, 
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are identically same, for consistency we must have 
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to the lowest order in r . In order to find out the physically relevant solution from 
the above two solutions in (1851) we argue in the following way: The functions Ri, R2 
etc are proportional to respective noise correlators (ITOl) in the coupled-KPZ equations 
( l69l) . Further, in terms of the original field variables u and b or h and 4>, if there is no 
cross-correlations, i.e., for To = the amplitude-ratio of the autocorrelation functions 
of u and b (or h and 0), A is unity. Since A = (1 — 7)/(l + 7) we then have 7 = 
when A — 1, i.e., in the absence of any cross-correlations. Thus we pick up that relation 
between 7 and T which goes to zero in the limit T goes to zero. Thus we write, 

i-r„ 



7 



.4 



2Tn 



i + r 

to the lowest order in 0(T Q ). Further, from its definition, r = (R x /Ri) 2 = N s in the 
lowest order. Hence, we obtain A = 1 — 2N S . Therefore, the relation (1HB1 agrees with 
what we find before from our DRG or SCMC calculations. 

The scaling exponents in the present coupled chain problem is identical to the 
single-DP problem; this is due to the identical nature of the functional flow equations of 
Ri with the corresponding single DP problem. Therefore, we obtain ( = 6/(8+ d) = 1/z 
(see also [17]) as obtained in our SCMC calculations before. 

As before in the DRG analyses of the problem, to complete our analysis here, it 
is required to demonstrate that the parameter r is marginal in the scaling limit and 
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can have arbitrary values. We begin by considering the flow equation for R x (x) at the 
RG fixed point. Noting that due to the odd parity nature R x (x) does not receive any 
fluctuation corrections we write 

^ = -£(3 - 40*,. (87, 

This yields, near the fixed point, 

Here, Cq is a constant of integration which is the value of R x (x) at small scale. The 
function i? x ( x ) is odd under x — > — x and is non-analytic at x = 0. As a result, within 
perturbative calculations, it does not receive any fluctuation corrections in the long 
wavelength limit. Hence even in that limit the value of R x (x) depends upon Cq, its value 
at the small scale. In contrast, the values of R\(x) and R2{x) in the hydrodynamic limit 
are independent of their values at small scales, since fluctuation corrections dominate 
over their bare values at large spatial scales. Therefore, the ratio T = [R x (x)/Ri(x)] 2 
at the large spatial scales, i.e., in the scaling regime, depends on Cq. The constant Cq 
has no fixed magnitude; it depends upon realizations of the disorder at small scales 
and hence can have arbitrary values. Therefore, the ratio [R x (x)/Ri(x)] 2 also can have 
arbitrary values in the hydrodynamic limit. This completes our FRG analysis. Our FRG 
approach to the problem, therefore, yields scaling exponents x = ~ § — f i - z = t + ff° r 
d — 1 and in the strong coupling phase at d > 2. It further yields the amplitude-ratio 
D b /D u = 1 — 2N S for d — 1 and in the strong coupling phase at d > 2. These results 
are in agreement with those from DRG and SCMC approaches above. 



5. Numerical analysis: direct approaches and lattice models 

5.1. Direct Numerical Solutions (DNS) of the model equations 

Having obtained several new results by the applications of three different analytical 
perturbative approaches on our model we now resort to numerical methods to 
supplement our understanding of the underlying physics from the above analytical 
approaches. In particular, we numerically solve (hereafter referred to as DNS) the 
model Eqs. and ( fTUl) in one and two dimensions by using pseudo-spectral methods 
with the Adams-Bashforth time evolution scheme [see Appendix (19.31) ]. Here, we 
consider only symmetric noise cross-correlations. We elucidate the scaling properties 
of the following equal-time correlation functions of /i(k, t) and 0(k, t): Chhik,0) = 
(h(k, t)h(— k, t)), C^k^O) = (0(k, t)0(— k, t)) in d = 1, 2. Since, as discussed before, 
the scaling exponents of the fields h and <fi are identical to each other, the ratio of 
the equal-time correlation functions A = C^(k, 0)/C , / 1 / l (k, 0) is a dimensionless number 
which is nothing but the amplitude ratio defined above. We examine the dependence of 
A on the parameter N s . We further consider the time dependence of the widths Wh(t) 
and W<f,(t) as defined above. In the statistical steady Whit) 2 and W^t) 2 approach 
the equal time steady state correlation functions C^(x = 0, 0) and C*00(x = O,O). 
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Therefore, Wh{t) /W^t) = y/A in the large time limit (i.e., in the statistical steady 
state). 

Before presenting our numerical results below we discuss a technical matter, namely, 
the measurement of the parameter N s . In our analytical work above, the parameter 
N s = (D s /D u ) 2 involves the ratios of the amplitudes of the cross-correlation function 
and the velocity auto- correlation function. Therefore, corresponding numerical works 
require measurements of the cross- correlation function as well, in addition to measuring 
the auto-correlation functions of u and b. It, however, is much more difficult to obtain 
data with sufficient statistics for the cross-correlation function amplitude since it is not 
positive definite. In view of this difficulty, instead of measuring N s we use its bare value 
as obtained from the amplitude of the noise cross-correlations in most of our analyses 
below. We denote this by N®. We would like to mention that the comparison of our 
numerical data with the already obtained analytical results will be largely qualitative, 
due to the reasons mentioned above. For our DNS studies the noises B\ and 9 2 in 
Eqs. ( fi~4"l) and (115jl are chosen to have correlations of the form 



such that D s (k) = -D s (-k) and D s (k) = D° s for k > 0. The resultant noise correlation 
matrix has eigenvalues D° ± D®. The fact that the noise correlation matrix should be 
positive semi-definite ensures that the upper limit of N® = (D®/D ) 2 is unity. Although 
N s depends monotonically on N®, due to the highly nonlinear nature of the equations of 
motion the dependence is not linear. We are able to measure N s only in lrfDNS below. 
Those measurements indeed show the monotonic dependences of N s on N®. We present 
our results in details below. 

5.1.1. Results in one dimension: In this section we present our numerical results from 
the DNS of the continuum model Eqs. ffT4"j) and (fT5j) together with the noise variances 
(I89p in Id. We have already found, from our analytical studies above, that the model 
Eqs. ([9]) and (flQj) together with the noise variances (jTTl) . (|T8|) and (TT9T) in Id yields 
scaling exponents x — ~ 1/2, i-e., Xh = 1/2 and z = 3/2. In addition, the amplitude- 
ratio Db/D u decreases monotonically with N s . Our results here confirm our analytical 
results as we describe below. 

We perform pseudo-spectral simulations of the model Equations in Id to examine 
the scaling behavior of the equal time correlation functions Chh{k,0), C^kyO), where 
A; is a Fourier wavevector, in the statistical steady states. The system sizes L chosen are 
L = 4096, L = 2048 and L = 6144 where L is the number of points in a one-dimensional 
lattice in real space. We present a log-log plot of the correlation functions versus k in 
Fig. [T] for L = 6144 (left) and L = 4096 (right); results from our runs with L = 2048 
have similar behavior. In all the plots the red point corresponds to the correlation 
function Chh and the green points correspond to the correlation function C^. The blue 
line with slope of -2 provides a guide to the eye for scaling regime in the plots with slope 



(^(k,t)^(-k,0)) = <0 2 (k,t)0 2 (-k,O)> = 2D 
(e 1 (k,t)^(-k,0)) = 2z J D s (k)5(t), 



5(t), 



(89) 
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log k log k 

Figure 1. Log-log plots of Chh(k) (red) and C r j )< p(k) (green) versus k for L = 6144 
with iV s °=0.9 (N s = 0.82) (top left) and N s = N° = 0.0 (bottom left), and for L = 4096 
with iV s °=0.9 (N„ = 0.81) (top right) and N s = N ( J = 0.0 (bottom right). The blue 
line, drawn as a guide to the eye, indicates a slope of -2, corresponding to the value 
Xh = 1/2 for the roughness exponent for the fields h and </>. Note that at finite N s 
there are amplitude differences between the scaling regimes of Chh{k) and C^(k) (top 
plots) where as for N® = 0.0 the amplitude differences disappear (bottom plots) (see 
text). 

—2 which corresponds to the roughness exponents Xh f° r the fields h and <p being 1/2. 
These values are exact results in the absence of cross-correlations and obtained in our 
one-loop DRG and SCMC above (Id) for finite cross-correlations. Our numerical results 
clearly yield a value of the roughness exponent Xh which is very close to the analytically 
calculated value. For our lrfDNS studies we estimate the parameter N s defined above by 
calculating the equal-time cross-correlation function for Fourier modes k in the scaling 
regimes and taking its ratio with Chh(k,0). In Fig. (CQ), for a given system size L, the 
amplitude differences between the scaling regimes of the correlators Chh and C^, which 
is same as the parameter A in Section 14.21 increases monotonically with N® (or with 
N s ), a feature which is in qualitative agreements with our analytical results above. 
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We now show the time-dependence of the widths Whit) and W^it) in Fig. ([2]) for 
system size L = 4096. We notice that the saturated amplitude difference between Whit) 
and W<j,{t) increases as increases from 0.0 to 0.9. Since the ratio of the saturated 
amplitudes of W h {t) and W^(t) yields the ratio A (W^/Wh = VA), we find that A 
decreases as N® increases. This is in accordance with the results as presented in Fig. 
dl]). Figure ([3]) shows the dependence of the parameter A obtained from the plots 




Figure 2. A plot of the widths Wh{t) (red) and W$(t) (green) versus time t for 
system size L = 4096 with iV s ° = 0.0, 0.9 from our Id DNS studies. Clearly the 
saturated amplitude difference increases as N® increases. 



above on N® for system sizes L = 2048,4096,6144. We find that the amplitude-ratio 
A decreases monotonically as N® increases in agreement with the behavior of our mode 
coupling analyses. We do not observe any noticeable systematic dependence of this 
behavior on system sizes. Therefore, results from our lrfDNS studies agree with those 
from the analytical studies mentioned above. 

5.1.2. Results in two dimensions: In this Section we present our results from our 2d 
DNS of the model Eqs. ( 1T41) and ( Tl5l) and compare with the analytical results already 
obtained. In our DNS studies (using pseudo-spectral methods) in 2d the system sizes we 
work with are 96 2 , 128 2 and 160 2 . Unlike in Id, the system exhibits a non- equilibrium 
phase transition from a smooth phase, characterized by logarithmic roughness (xh = 0) 
and a independent of iV°, to a rough phase characterized by algebraic roughness (xh — 0) 
and a decreasing A as N® increases. In the smooth phase, A is fully determined by the 
bare ratio D®/D®. Therefore, a simple way of ascertaining which phase the system is 
in, is by measuring A. 

In our 2c£DNS studies there two tuning parameters to reach the rough phase of 
the system. They are D® and iV°. The crossover from the smooth-to- rough phases 
is formally determined by the value of the coupling constant G = A-D 2 /V 3 (see above). 
Here the renormalized parameter D u has a monotonic dependence on the bare amplitude 
D®. For small D®, renomalized G is zero and the system is in its smooth phase. With 
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Figure 3. A plot of the parameter A versus N® for our Id DNS studies for system 
sizes L = 2048, 4096, 6144. The general trend that A decreases with an increase in N® 
is observed for the all system sizes we worked with. 

increasing D® the system eventually crosses over to the rough phase where various values 
of a have been obtained by tuning N®. For each system size we obtain the widths Whit) 
and Wtfrif) as a function of time t. As for the height field in the well-known KPZ equation 
plots of Whif), as a function of time, have a growing part and a saturated part. The 
system size dependences of the saturation values of the widths Whit) and W^t) yield the 
values of the roughness exponents of the corresponding fields. Since we are interested in 
the statistical properties of the rough phase, in our 2<iDNS studies we access this phase 
by sufficiently large D® for reasons as explained above. As shown below, we find A < 1 
for non-zero N®. This ensures that we are indeed able to access the rough phase in our 
DNS runs. We present our results from system sizes 96 2 and 128 2 in Fig. (j3J) for two 
values of the parameter N® = 0.0, 0.9 for each system size. The data from 160 2 runs 
show similar behavior. We determine A and the ratio Xh/ Z - Note that for each system 
size the amplitude differences differences between the saturation values of the widths 
Whit) and W^it) (i.e., the parameter A in Sections 14.21 and 14.31) increase monotonically 
with N®, in agreement with our SCMC and FRG analyses. Such nontrivial dependence 
of A on Ng is a key signature of the rough phase. From our DNS studies with system 
sizes 96 2 , 128 2 , 160 2 we calculate the parameter A. We show its dependence on N® in 
the Fig. (J5j). As in Id we see that A decreases monotonically as N® increases from zero. 
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Figure 4. Plots of Whit) (red) and W${t) (green) versus t in the rough phase for 
the system sizes 96 2 (top) and 128 2 (bottom) for iV° = 0.0, 0.9. Note the increasing 
differences between the saturation values of Wh and W^ with increasing . 

We now determine the scaling exponents Xh and z. Due to the very time consuming 
nature of the 2d simulations the quality of our data for 2d are much poorer than 
those obtained from our Id runs. Because of this difficulty we do not extract the 
roughness exponent directly by plotting different correlation functions in the steady 
state as functions of wavevector k. Instead we obtain the ratio Xh/ Z by plotting the 
width Wh(t) as a function of time t in a log-log plot for t < saturation time scale. The 
slope yields the ratio Xh/ Z ', we obtain Xh/ Z — 0.2 ± 0.1 which is to be compared with 
our analytical result 1/5, obtained by means of SCMC and FRG calculations. At d = 2 
SCMC yields z = 1.66 which is close to Ref. [29] whereas our DNS yields z ~ 1.60 ±0.1. 
Numerical studies of our type, performed on bigger system sizes, should be able to yield 
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Figure 5. A plot of A versus N® for various system sizes from our 2dDNS studies. 

highly accurate values for the scaling exponents which could be compared systematically 
with our SCMC/FRG results and those of Ref. [29]. However, we refrain ourselves from 
making such detailed comparisons due to the rather small system sizes we have worked 
with. We present our result in Figj6] below for the system sizes L = 128 2 and L = 160 2 
with N® = 0.0 and N® = 0.9 respectively. Note that, within the accuracy of our 
numerical solutions, the ratio Xh/z does not depend upon the value of N®, suggesting 
that the scaling exponents are independent of symmetric cross-correlations in agreement 
with our SCMC (Section H3D and FRG (Section Oj) above. 

We conclude this Section by summarizing our results. We performed DNS of the 
model Eqs. ( fT4l) and ( fl5l) at d = 1, 2 in the presence of the symmetric cross- correlations 
parametrised by N®. By tuning the bare amplitude D we are able to access the rough 
phase in 2d. In both dimensions we find that in the rough phase the ratio D^jD u 
decreases monotonically from unity as N® increases from zero. In Id we find Xh — 0.5 
very accurately which is in good agreement with the analytically obtained results. In 
2c? we find Xh/z = 0.2 ± 0.1 which is to be compared with our analytical estimate of 
1/5. Extensive DNS studies on larger system sizes would be required to calculate the 
scaling exponents with high accuracy. 
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Figure 6. Log-log plots of Chh{k, t = 0) and C^{k,t = 0) versus k from our 2dDNS 
studies: (i) L = 128 2 , iV s ° = 0.0, (ii) L = 160 2 , N° = 0.9. The slope - 0.2 yields Xh/z 
(see text). 



5.2. Monte Carlo simulations of the coupled lattice-gas models in one dimension 

So far in the above we have studied the applications of several methods, analytical as 
well numerical, on the continuum model Equations and (fTU|) [or Eqs. (fl4|) and ([15]) ] 
and obtained several results concerning universal properties of the statistical steady- 
state of the model. Note that all techniques used above have been applied on the same 
model. To complement our studies we discuss the results from the lattice-gas model in 
this section which we constructed [see Section [3], and compare with our results already 
obtained above. 

In this section we simulate our proposed one-dimensional lattice-gas models for 
the model Eqs. fl2Tj) . By using relations fl23|) we calculate the ratios of appropriate 
correlation functions in the steady state. We use Monte Carlo methods to simulate 
our models. We extended the Restricted-Solid-On-Solid (RSOS) algorithm [20] and the 
Newman-Bray (NB) algorithm [22] for the KPZ surface growth phenomena to construct 
the coupled lattice-gas models in Id for Eqs. (l2Tj) with cross-correlations. In such 
models particles are deposited randomly from above and settle on the already deposited 
layer following certain growth rules which define the models. One typically measures 
the widths of the height fluctuations of the growing surfaces as functions of time t. Note 
that our numerical works on the lattice-gas models are restricted to lattices of modest 
sizes (L = 4096 is the largest system size considered). This is due to the fact that 
the numerical generation of noise cross-correlations of the type we considered is not 
uncorrelated in space in Id, rather it has a variance proportional to 1/x. Generation 
of such noises is a very time consuming process: we first generate the noises in the real 
space. These are uncorrelated Gaussian random noises without any cross-correlations. 
We then bring them to Fourier space by Fourier transform. Inverse Fourier transforms 
of particular linear combinations of them yield noises with finite cross-correlations which 
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we use in our studies. Although we have used Fast Fourier Transforms, they are still 
computationally rather time consuming and hence reduce the over all speed of the code. 
In contrast, the more common lattice-gas studies on the KPZ equation do not require 
Fourier Transforms of the noise (since they do not have any noise cross-correlations) 
and hence are much faster. This allows them to study up to much larger system sizes. 



5.2.1. A coupled lattice-gas model in one dimension with RSOS update rules: The 
restricted solid-on-solid (RSOS) update rules involve selecting a site on the lattice 
randomly and permitting growth by letting the height of the interface at the chosen 
site increase by unity such that the height difference between the selected site and the 
neighbouring sites does not become more than unity [20]. Our model involves two 
sublattices where the height fields on each of the sublattices satisfy each of Eqs. (121)) . 
In our coupled lattice model each sublattice is evolved according to the RSOS update 
rule described above. The selection of the random sites in the two sublattices may be 
correlated. This correlation models the noise cross-correlations of the model continuum 
Eqs. © and ([lO]) or (EH). 

We simulated our coupled lattice-gas model described above and calculate the 
widths Whif), W<p(t) of the growing surfaces h(x,t) and 4>(x,t) as functions of time 
t. Similar to the single- component KPZ equation, the plots have two distinct parts - 
an initial growing part and a late time saturated part. Below we present our results 
graphically: We plot Whit) (red) and W^{t) (green) versus t for system size 1024 in 
Fig. ©. Results from system sizes L = 2048, 4096 show similar behavior without any 
systematic system size dependence. 
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Figure 7. Plots of Wu{t) (red) and W^m (green) versus t for system size L = 1024 
for iV° = 0.0, 0.9 in our coupled lattice-gas model based on the RSOS algorithm. Note 
the increasing difference between Wh and W<j, with increasing N°. 



In Fig. [H]we present a plot of A versus N® for the system sizes L = 1024, 2048, 4096. 
As in our results obtained analytically and DNS, A decreases monotonically from unity 
as N® increases from zero. We do not observe any systematic system size dependence. 
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Figure 8. A plot of A versus N® for various system sizes obtained from numerical 
simulations of our coupled lattice-gas models with the RSOS algorithm. 

z to complete the discussions on the universal properties of the lattice-gas model. Below 
we present a log-log plot of the equal-time correlation function Chh(k) = (h(k, t)h(—k, t)) 
in the steady state as a function of the Fourier vector k [Fig. We find that slope in 
the scaling regime (small k) is very close to -2 corresponding to Xh — 1/2 as predicted 
by analytical means before. 

An alternative way to obtain the scaling exponent Xh is by finding the dependence 
of the saturated value of the widths Whit) and W^{t) on the system size L. Since, after 
saturations in the steady states Wh,$ ~ L 2xh with L being the system size, we obtain 
Xh from the slope of the plots of logarithm of the saturated values of the widths versus 
logarithm of the corresponding system sizes. We find Xh = 0.49±0.05 which is very close 
to the analytically obtained value, Further, from the time-dependences of the widths 
Wh{t) and W^it) the exponent-ratio Xh/z can be obtained. Fig. HU1 shows a log- log plot 
of Whit) and W^it) versus t from a lattice size L = 16384. For such a large lattice 
size the saturation time is very large and hence we show only the growing part of the 
curves. The slope yields the ratio Xh/%- The red and green points refer to Whit) and 
W^it), and the blue line a slope of 0.32 indicating Xh/z = 0.32 which is very close to the 
analytically obtained value of 1/3. In short, therefore, the results from the Monte-Carlo 
simulations of our coupled lattice-gas model in lei, based on the RSOS algorithm, yield 
results in close agreement with those obtained through analytical means and lciDNS. 
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Figure 9. A log-log plot of the correlation function Chh{k) versus k for system size 
L = 2048 from our coupled lattice-gas model with the RSOS algorithm. The straight 
line with a slope of -2, corresponding to Xh = 1/2, represents the scaling regime of 
Chh(k) for small k (long wavelength limit). 



5.2.2. A coupled lattice-gas model in one dimension with the NB update rule: As 
mentioned before, this model uses the mapping between the KPZ surface growth problem 
and the equilibrium problem of a directed polymer in a random medium (DPRM) which 
are connected by the non-linear Cole-Hopf transformation /i 12 (x, t) = (2^/A)lnZ lj2 
leading to the Eqs.l fTTT) for the partition functions for the two DPs. Further, one 
uses the following update rules for h\$. 

kit) = h l (t) + A 1 ^Q(t), 

hi{t + A) = kit) + (v/X) In [l + (Au/a 2 )i: jnni [e x ^-' h ^ - 1]] , (90) 

where A and a are the grid scales for time and space, respectively, and jnni indicates 
sites j and i being nearest neighbors. Then, taking the strong coupling limit A — > oo 
one finally obtains 

hi(t + A) = max jnni (hi, {hj}). (91) 

This is same as the zero-temperature DPRM algorithm, written in terms of the fields 
ft<i,2- We implement the above growth rule in Id with system sizes L = 1024, 2048, 4096. 
Functions Whit) and W^(t), as defined in Sec J5.2.11 are obtained by taking appropriate 
linear combinations of the equal-time correlators of h 12 . We present our results 
graphically in Fig. (II ip for system size L = 1024 below. 
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Figure 11. Plots of Wh(t) (red) and W$(t) versus t for system size L = 1024 for 
various N s in our coupled lattice-gas model based on the NB algorithm (see text). 
Note the increasing difference between W h and with iV s ° = 0.0, 0.55, 0.9. 



From Fig. (fTTj) it is clear that after saturations, the differences between Wh and 
Wtj,, measured by A, increase with N®. Below in Fig. (TI2]) we show a plot depicting 
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the variation of A versus N® obtained from the simulations of our coupled lattice-gas 
model with the NB update rule in Id. As before, we find that A decreases monotonically 
from unity as N® increases from zero. We further determine the scaling exponents Xh 
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Figure 12. A plot of A versus N® obtained from numerical solutions of our coupled 
lattice- gas model with the NB update rule for system sizes L = 1024, 2048, 4096 in 
Id. 

and z from our numerical results, as we did above in Section 15.2.11 Fig. (1131) shows a 
log- log plot of Whit) and W$(t) versus t for the system size L = 4096. We find that the 
initial slope is 0.32 ± 0.02 which is very close to the analytically obtained Xh/z — 1/3. 
Furthermore, the equal-time correlation function Chh{k, 0) as a function of k is shown in 
a log- log plot [Fig. [H]. We find that in the scaling regime (small k) k~ 2 is an excellent 
fit for Chh{k) yielding x — 1/2- This, together with the fact that Xh/ Z — 1/3 give 
z ~ 3/2, close to what we estimated by other means in Id as reported above. 

We close this section by summarizing our results from the Monte-Carlo simulations 
of the two lattice-gas models proposed by us. Both the models, based on the RSOS 
and the NB algorithms respectively, exhibit similar dependence of the parameter A over 
Ng. They further yield the same scaling exponents. All these are in close agreements 
with our analytical results obtained from DRG/SCMC/FRG methods as well as lrfDNS 
studies. 



Scaling and universality in coupled driven diffusive models 



10 




log k log k 



Figure 13. Log-log plots of Wh(t) (red) and W${t) (green) versus t for system size 
L = 4096 with N° = 0.99 (left) and N° = 0.0 (right) respectively. These are obtained 
from our lattice-gas model with the NB algorithm The initial slopes are very close to 
the analytically obtained value (sec text). 
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Figure 14. A log-log plot of Chhik) versus k from simulations of our coupled lattice- 
gas model with the NB algorithm (L=4096). The black line is a guide to slope -2 which 
agrees very well with our data in the scaling (small k) region. 



Scaling and universality in coupled driven diffusive models 



41 



5.3. Lattice-gas models and the continuum equations of motion: Universality for 
systems out of equilibrium 

In the above we have obtained results on the scaling exponents and the correlation 
functions in the statistical steady states from different approaches. In particular we 
compare our results in Id on the dependence of the amplitude ratio A on the symmetric 
cross-correlations N® from the numerical solutions (DNS) of the continuum model Eqs. 
([9]) and ( flQl) . and the numerical simulations of the lattice-gas models based on the RSOS 
and NB algorithms. Since in none of the cases we considered we find any systematic 
system-size dependence, we take average of our results from different system sizes for the 
same value of iV° for a given type of study (DNS/NB/RSOS). We present our results in 
a composite plot ffT5]) below. We find that the general trend that the ratio A decreases 




s 



Figure 15. A plot showing the variation of the amplitude ratio A as a function of the 
symmetric cross-correlations N® (see text) in Id from RSOS (red x), NS (green *) and 
DNS (blue diamond). 

from unity when N® increases from zero persists in all the cases. Moreover, all the cases 
yield values of the scaling exponents which are very close to each other and in close 
agreement with the analytically obtained values as well. 
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6. Summary and outlook 

We have considered the statistical properties of the generalized Burgers model (GBM) 
egiven by equations (Q and fflQj) . driven by conserved zero- mean Gaussian noise sources 
with variances ( TPTl [TBI) and ( [19]) . We have undertaken a variety of analytical and 
numerical investigations to uncover the universal properties of the correlation functions 
in d dimensions. 

Upon employing analytical approaches like DRG, FRG and SCMC we have 
calculated the scaling exponents and the amplitude ratios at d — 1 and at d > 2. Each 
of these methods have a certain range of applicability and also limitations. The DRG 
method yields results only for the strong coupling phase in d = 1 and the unstable critical 
point (for d > 2), for both symmetric and antisymmetric cross-correlations. The FRG 
method works for the strong coupling phases at d — 1. For d > 2 it is able to capture 
the physics in the strong coupling phase, but only for symmetric cross-correlations. The 
SCMC method is useful to study the strong coupling phases at d = 1 and d > 2 for 
both symmetric and anti-symmetric cross-correlations. Taking all the methods from the 
various analytical methods together a coherent picture emerges. We broadly find that 
symmetric cross-correlations affect the amplitude ratio and the anti-symmetric cross- 
correlations affect the scaling exponents. Since, in general, both symmetric and the 
anti-symmetric cross-correlations are likely to be present, one may encounter scaling 
exponents and amplitude ratios varying with the strength of cross-correlations. In 
addition to our analytical studies, we have performed extensive numerical work in Id. 
We have numerically solved the continuum model equations using the pseudo-spectral 
method in one and two dimensions, and numerically simulated two equivlent lattice- 
gas models in Id. In all the cases we numerically calculated the equal-time correlation 
functions of the two height fields h and 0, and measured their spatial scaling exponents 
and the amplitude ratios. Our numerical results on the scaling exponents agree very well 
with those obtained analytically for d — 1. Further, our results on the dependence of 
the amplitude-ratio on the strength of the cross-correlations from the numerical studies 
agree broadly with those from the analytical studies. 

In 2c? we have numerically solved the continuum model equations (fl4j) and (IT51) 
in the presence of symmetric cross-correlations only. The amplitude ratios display 
dependences on the cross-correlations similar to those obtained analytically. We also 
obtained scaling exponents numerically. Due to the difficult nature of the 2d simulations 
these are less accurate than their Id counterparts. Since we use a pseudo-spectral code, 
our method is free from the controversial issue of real-space discretization [35J. Owing 
to the small system sizes and the high saturation time for the widths, the quality of our 
data are less satisfactory than those obtained from our Id models. We are only able 
to measure the ratio Xh/z from our numerical studies in 2d. Nevertheless, the effects 
of noise cross-correlations have qualitatively the same effects on the amplitude ratio. 
Detail studies to yield the scaling exponents in 2d with high accuracy will be presented 
elsewhere 1361. 
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Our results potentially have broad implications for natural systems with coupled 
degrees of freedom. There are many such examples - MHD turbulence, turbulence in 
a binary fluid mixture, biologically relevant systems, e.g., coupled TASEPs (Totally 
Asymmetric Exclusion Process) modeling transports in cellular systems, dynamics of 
motors and microtubules, dynamics of visco-elastic active polar gels etc. The general 
conclusions of our investigations, namely, the relevance of the strength and the symmetry 
of the cross-correlations in determining the long wavelength properties, are expected to 
be visible in such systems as well, although detailed accounts are likely to depend on the 
systems in question. To test our predictions experimental measurements should include 
correlation functions and effective transport coefficients. These are likely to be affected 
by the presence of cross-correlations. In driven fluid-like turbulent systems (MHD, 
binary fluid mixtures) cross-correlations can be experimentally controlled by tuning the 
external (stochastic) forces. In such systems one needs to measure correlation functions 
of velocity and magnetic fields (in MHD) or velocity and concentration gradient 
fields (in binary fluid mixtures) in the scaling limit and obtain the scaling exponents 
characterizing the correlation functions in the hydrodynamic limit as functions of the 
amplitudes of the cross correlation functions. 
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8. Glossary 

For the convenience of the reader we provide below a list of symbols for the functions 
and parameters which are used in the bulk of the article. 

• u(x, t), b(x, t) =>- Burgers fields. 

• = u ± b Elsasser variables. 

• u = V/i, b = V0. Fields h and are (non-conserved) height variables. 

• 7^ = V/ii,2, hi s 2 are (non-conserved) height fields; hi^ also represent the free 
energies of two (identical) directed polymer in a random medium; /i 12 = \{h ± (/)). 

• z/o, fio =r- Bare (unrenormalized) viscosities for the fields u and b; v, fi =>- 
Renormalized fluid and magnetic viscosities. 

• Ai, A2, A3 =>• coupling constants. They are set to unity in most of what follows. 
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(j U) &(k, uj) =>■ Renormalized propagators/response functions of the fields u and b. 

S(k, uj) => Self-energy function of u and b; Y =>- amplitude of £(k, 0), Y = u = ji 
in the present model. 

CJi, C^a, Cf: =>• Auto correlation functions of u and b, and cross correlation function 
respectively, i,j =>■ Cartesian indices. 

D®, D® =>■ Bare amplitudes of the auto-correlations of the external Gaussian noise 
sources in the u and b equations [Eqs. ( |9ll9l) ]: D u , Df, =>- corresponding renormalized 
(scale- dependent) amplitudes. 

D®, D® => Bare amplitudes of the symmetric and anti-symmetric parts of the noise 
cross-correlations in Eqs. ( |9f9l) ; D s , D a =>- corresponding renormalized amplitudes. 
P m = u/fi =>- Renormalized Prandtl number. 
A = Df,/D u =>- Dimensionless renormalized amplitude ratio. 

N° = (Dgf -D.°) 2 , N s = (D s /D u f, N a = [D a /D u f. 
D = D° U + D° b , D = D° U -Dl 

V\,V2 Random potentials embedding directed polymers (DP) A and B. 
Zi, Z2 =>- Partition functions of DP A and B. 

Ri , R2 , R, R Different variances of the random potentials V\ and V%. 

Jrj^y- Since, -Ri(x), i?(x), i? x ( x ) have same scaling behavior in 
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the hydrodynamic limit, 7, To are constants. To the lowest order To = N s , 7 = A. 

G = x2 Ji u =>- Dimensionless coupling constant in one-loop dynamic renormalization 
group calculations in the presence of symmetric cross-correlations; G a = A J^" =>- 
the same for anti-symmetric cross-correlations. 

T =^> Vertex generating functional. 

e = 2 — d where d is the physical space dimension. 

Xui Xb, z =>- Two roughness and dynamic scaling exponents respectively of the fields 

u and b. In the model here, Xu = Xb — Xi X + z = 1- 

Xhj X<t> =^ Roughness exponents of h and (f>; Xh = X4>-, Xh + z = 2. 

( = 1/z => scaling exponent characterizing the scaling of transverse fluctuations of 

directed polymer with its longitudinal length. 

T — 1/z/q =>- Temperature of a DP in a random medium. 



9. Appendix 

9.1. Symmetries and Ward identities 

In this section we elucidate the continuous symmetries under which the equations of 
motion and ffTOl) remain invariant. These allow us to construct exact relations 
between different vertex functions which in turn impose strict conditions on the 
renormalization of different parameters in the model. 
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We begin by rewriting the Eqs.([9"l) and fllOp by expressing the Burgers b- fields 
(whose mean value is always kept zero here and below) as a sum of a constant vector 
and a space-time dependent part, i.e., splitting the fields b(x, t) — > b(x, t) + B0/A2 
where B0/A2 is a constant vector: 

+ VB ■ b + Ax Vu 2 + A 2 V6 2 = u V 2 u + f (92) 

at 

and 

db 

+ (A 3 /A 2 )VB ■ u + A 3 V(u • b) = // V 2 b + g. (93) 

In this notation the mean Burgers b-field is (b) + B = 0. 

The generating functional corresponding to Eqs. ([9]) and (fTUj) is written as [37J [38] 
jij J2, j2] = / DuiDuiDbiDbi exp[<S] [39], where the action functional S is given as 



S[$\ = - J bibj D ij- J UiUjD^-i J u i b j D Xs -i J u i b j e ijp k p iD Xa 

- % J Ui[d t Ui + y WU 2 + y Wb 2 + VB • b - l^jVV] 

- % [ h[dA + A 3 V(u ■ b) + (A 3 /A 2 )VB • u - m V%]. (94) 



From the dynamic generating functional one obtains the vertex generating 
functional T [39J through 

fKfi„M ( ] = W^hlJi - J - JjA - //A -j Ik, (95) 

where W = In Z is the generator of the connected diagrams and T is the generator of 
the one-particle irreducible diagrams [40J. 

The equations of motion (Q and (fTUl) or the action functional given by (j94|) are 
invariant under the following transformations: 

• The Galilean transformation (TI): u(x, t) — > u(x + Uot, t) + Uo, §1 ~^ §1 ~ Auo-V, 
and b — >• b [HJ [15l HTJ with Ai = A 3 = A in Eqs. ( |92|) and (|93l This implies, as we 
show below, non-renormalization of Ai [HI [TJ3, [23] . 

• The transformation Til: Bo — > Bo + A 2 5, b(x, t) — > b(x, i) — 5, u — >• u. Here the 
shift S is a vector. 

Let us first consider the invariance under the Galilean transformation TI and the 
corresponding Ward identities. The invariance of the vertex generating functional T 
yields 

ST = u - / dti\qt[———u j (q,t) + —- — -Uj(q,t) 

+ 4^ Mq ' t) + w' J(q ' t,, 



u - / [ dtS(q) r = 0. 



(96) 
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By taking functional partial derivatives with respect to different fields we obtain (in 
time space) 

K<tf + <?i'*")f u jUk = - J dtt ajUiUk (q', -q', t"; q, t) | q=0 , (97) 

K<tf + q"t")Tu 3 u k = ~ J dtf aj u kUi (q', t'; -q', q, t) | q=0 , (98) 

Kit' + W)h ]hk = -ydA. Mi (q',t';-q',t";qt)| q=0 , (99) 

A(^ + gfOf^ = "/ ^^(q',t';-q',t";q,t)| q=0 . (100) 

The above Ward identities may, equivalently, be written in the frequency space as 
d ~ 

A%^r fiiUfc (q,o;;-q, -u) = r fiiU . Ufc (q, u; -q, -w; 0, 0), (101) 

Ag^r fiiiifc (q,a;;-q, -w) = r fi . fifeU .(q, w; -q, -w; 0, 0), (102) 
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-a;) 


— r~ 
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— UJ 
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= f - - 
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— UJ 


;0,o), 


-q 


- W ) 


— f- ( 

bjb k Ui y 




-q, 


-u; 


0,0), 


-q, 


-a;) 


— f- - { 

bjb k Ui 1 




-q, 




0,0). 



<9 ~ 

^ftg^r Sj . 6fc (q,a;;-q,-w) = r g . 6tU4 (q, cj; -q, -w; 0, 0), (103) 
d ~ 

A %^r^^(q,^;-q,-^) = r^ SiU .(q, w; -q, -w; 0, 0). (104) 



Further, as discussed in Ref. [22], t ne diffusive dynamics of the model and the 
corresponding q dependences of the vertices yield the exact relation 

f« i% .(q = 0,w) = zw, f Si6i (q = 0,a;) = iw. (105) 

To proceed further we first define the renormalization Z-factor for a quantity k in the 
usual way: Z k kr = k where the suffix R refers to renormalized n. These then lead to 

Z u Z a = l, Z b Z~ b = l. (106) 

For the choice Z u = 1 (one of the renormalization Z-factors can be set to unity without 
any loss of generality) we then obtain from the above Ward identities that the non- 
linearity A does not renormalize. 

We now discuss the symmetry TIL This symmetry has recently been discussed in 
the context of a symmetric binary fluid mixture model [12]; see also [13] for a similar 
symmetry in a different physical context. The essence of the invariance under the 
transformation Til is the following: We decompose the total Burgers b-field as a sum 
of a constant vector B and a space-time dependent part b(x, t) while writing down 
the Eqs. (j92|) and fl93]) where B is completely arbitrary. It should be noted that B is 
not the mean b-field, i.e., Bo is not the average of the total b-fields. The latter, in our 
notations, is given by B /A 2 + (b(x, £)) since B /A 2 + b is the total field. In fact, in 
our calculations here and below we have set the mean magnetic field to zero which can 
be ensured by adding appropriate counter terms in the action above. We did not show 
them explicitly, but these are built in our calculations below. 

Under the transformation Til the mean Burgers b-field B0/A2 + (b(x, t)) remains 
unaffected. We argue that such a freedom to decompose the total magnetic fields into 
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two parts should hold for the renormalized versions of Eqs. (1921) and (|9"5j) as well. In other 
words, in the renormalized versions of Eqs. fl92|) and fl93|) one would be able to combine 
the constant and space-time dependent parts of the Burgers-b field and express them in 
terms of the total Burgers-b field. The resultant renormalized equations of motion must 
be identical with those which are obtained from the bare Eqs. ( 1921) and ( 1931) by first 
writing them in terms of the total Burgers-b fields (by combining B and b) and then 
expressing them in terms of the renormalized fields and parameters. This can happen 
only when the second and the fifth terms of Eq.( l92i) renormalize in the same way. We 
arrive at the same conclusion below by using a more formal language. Before closing 
this discussion we again emphasize that here we consider a homogeneous and isotropic 
system, i.e., there is no mean Burgers-b field in the system. Furthermore, to clarify the 
matter from a technical point of view we note that the imposition of the condition of 
zero mean Burgers-b field can be achieved by any of the two choices below: 

• The choice (b(x, t)) = —Bo ^ 0, such that the actual value of the mean Burgers-b 
field (b(x,t)) + B = 0, 

• The choice (b(x, t)) = = B . This too, of course, maintains zero value for the 
mean Burgers-b field. 

Note that the choices above are connected to each other by the transformation Til of our 
manuscript. We now demand the physical requirement that the renormalized equations 
of motion (and hence all measurable quantities) would be the same, regardless of the 
choice that one may make to obtain them. Such a requirement is nothing but the 
invariance under the transformation TIL In particular, we work with the second choice 
above, namely, (b(x, t)) = = B in our perturbative calculations, which is the most 
convenient one. However, the fact that our perturbative RG is in agreement with the 
Ward identity originated from the transformation Til, ensures that the results from our 
perturbative RG are actually independent of the choices mentioned above. 

To proceed further, we note that in the present problem one has Z u = = 1 
[23| |4*4"] . Since the corresponding vertex generating functional is also invariant under 
the transformation Til we find 

X2 Jb^ = 5bi(q = o,u = oy (107) 

This leads to 

2 SB 0i 5b j {k)Su j (-k) ~ 56 i (k)5u J -(-k)56 i (q = 0)' ^ ' 

Therefore, for the renormalized action to be invariant under the same transformations, 
we must have 

Z Bo Z b = Z X2 Z 2 h . (109) 

We assume the vector Bo to have components Bq^ where /i is arbitrary. From the 
action functional (1941) then, 
5 5 2 T 



5B 0fl Sb^Suii-k) 



ik^Sij. (110) 
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Therefore, its renormalization Z-factor must be unity leading to Zb Zi, = 1. Due to the 
Ward identity (11081) the two sides of it renormalize the same way and the Z-factor of 
the right hand side is Z^ 2 (we have used the fact that Z„ = 1). The left hand side has 
an overall Z-factor given by Z\ 2 Z B ~^Z b ~ 1 = Z\ 2 where we have used the relation (I109[) . 
This yields, 

Zx 2 Z 2 b = 1. (Ill) 

We, therefore, conclude that the coefficient of the effective vertex constructed by using 
V^b does not renormalize. Hence, A2 can be set to unity by treating all the Burgers-b 
fields as the effective Burgers-b fields "v/Xjb. In fact, non-renormalization of A2 can be 
argued in a less formal way: Let us assume that under mode eliminations and rescaling 
A2 — > A2/A. This scale factor of 1/A can now be absorbed by redefining the units of the 
Burgers-b fields by b — > v^/Vb. Mathematically this would mean Z\ 2 Z 2 = 1 - the same 
conclusion as in Eq. flllip obtained by using a more formal way. Since the Eq. ( fTOl) is 
linear, redefining (rescaling by a multiplicative factor) b-fields leaves it unchanged (see, 
e.g., Ref.jH]). Henceforth, all b fields here are to be understood as effective fields with 
the coefficient of the effective vertex being set to unity. It should be noted that under 
the rescaling b — > v^Ab the stochastic force Q is also scaled by a factor y/X. However, 
this does not affect any of the conclusions as the assignment of canonical dimensions to 
various fields and parameters is done after absorbing A2 in the definition of b (or, by 
assigning zero canonical dimension to A2). 

9.2. Noise cross- correlations in the real space 

We now establish the form of the noise cross-correlations in real space in one dimension. 
Note that unlike the auto-correlation functions cross-correlations are not squares of 
modulii of functions; rather they are products of two different complex functions in the 
Fourier space. Hence they can be positive or negative, real or imaginary in the Fourier 
space. This, together with the fact that noise cross- correlations are real in the direct 
space, and the properties of the fields u and b under parity inversion yields that noise 
cross-correlations in Id [see Eq. (l26l) ] are imaginary and odd in Fourier wavevector k. In 
Id it has the form iDk/\k\. An inverse Fourier transform yields the form in real space 
D(x): 

OO 1 POO POO 

D(x) = I dke tkx ib— = if) / dke ikx - if) / dke~ ikx 



00 IN Jo Jo 

poo 

2D / sm(kx)dk. (112) 







The last integral above can be evaluated including a convergence factor (equivalently 
by using contour integration along a closed contour from the origin to 00 along the real 
axis, then along a circle at infinity anticlockwise up to an angle tt/2 and finally back to 
the origin along the imaginary axis). This yields 



POO 

D(x) = -2bLt^ i sm(kx)e-' yk dk (70 > 0) 
Jo 
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= - 2DL£h j _k)9 — - — = -2D-. (113) 

ix — 70 x 

Therefore, the cross-correlations have a form 1/x in real space in Id. This has the same 
dimension as the <5-function in Id, but has a longer range than the 5-function. 

9.3. Pseudo-spectral schemes for DNS 

We would like to solve the stochastically driven equations (Q and (fTOj) in a one- 
dimensional box. We use a pseudo-spectral scheme with periodic boundary conditions 
(PBC). Pseudo-spectral schemes have a long history in the numerical studies of fluid 
turbulence [see, e.g., Ref. [15]]. The PBC are the simplest boundary conditions ensuring 
that there are no surface effects on the bulk properties that we propose to measure. They 
ensure that the fields u and b can be expanded in terms of discrete Fourier modes, labeled 
by k. The highest fc-mode will determine the smallest scale that can be resolved in the 
simulations. In pseudo-spectral schemes spatial derivatives are evaluated in the Fourier 
space, while products of fields are evaluated in the real space. In terms of the Fourier 
modes u{k) and b{k) the equations (Q and (jTUl) take the form 

—+— t)u(k-q, t)+— t)b(k-q, t) = -uk 2 u{k, t)+f, (114) 

i 

db 

-^ + ik} j u(q, t)b(k - q,t) = -fik 2 b(k, t) + g. (115) 

Q 

We consider the case v = fi. For time integration we use an Adams-Bashforth scheme 
(step size St) 

exp(-2vk 2 5t) - ^J fc n - 1 exp(-2z/A; 2 5t) + /£ exp(-isk 2 5t), (116) 



and 



K+1 ~ b k ex v(-vk 2 6) _ Z E n eM _ uk 2 5t) _ \ E n k - l exv{-2vk 2 5t) 



St 2 '2 

+ g n k ex V {-vk 2 St). (117) 

Here, Q = l/2u 2 , J = 1 / 2£> 2 , E = ub which are algebraic products of functions in 
real space, superscripts n etc refer to the n-th time-step. Functions fl, J and E are 
evaluated in real space which are then brought to Fourier space by Fourier transforms. 
Note that, unlike in the commonly used Euler method of time integration, here one 
solves the linear parts of the Eqs. of motion exactly and one requires the solutions at 
steps n and n — 1 to solve at the (n + l)-th step. Time evolutions of u and b takes place in 
Fourier space. Fast Fourier Transforms (FFT) are used to perform Fourier transforms. 
We, in particular, used FFTW (http://www.fftw.org) routines for this purpose. 



Scaling and universality in coupled driven diffusive models 



50 



The choice of updating time step 5t:- For such nonlinear equations as ours there is no 
analytic method available to choose the updating time step 5t which ensures the stability 
of the scheme. To get a measure of stability criteria we perform a von-Neumann linear 
stability analysis on the linear parts of our PDEs in /c-space leading to the following 
condition on 5t: 

|1 - uk 2 5t\ > 1 2 > uk 2 5t. (118) 

For our run with L = 6144 for a linear box size 2n, k max = 3072. We choose 
v = 5 x 10~ 9 , St = 0.008. Thus vk 2 bt = 0.004 < 2. We start from random initial 
conditions on the u and the b fields, allow it to reach steady states (we monitor 
it by noting the time evolutions of the energies of the u and b fields, which should 
fluctuate about their mean values in the statistical steady state). We then perform all 
our measurements in the statistical steady state obtained. 

9.4- Generation of noises 

For our DNS studies in Id we need to generate stochastic forcings / and g with variances 
(in /c-space) 

(f(k,t)f(-k,0)) = 2Dk 2 5(t), 
(g(k,t)g(-k,0)) = 2Dk 2 5(t), 

(f(k,t)g(-k,0)) = 2Wk\k\S(t). (119) 

To obtain / and g we first generate two independent sets of random numbers with 
Gaussian distributions with zero mean and of variances 2(D + D)k 2 and 2{D — D)k 2 . 
Appropriate linear combinations of these random numbers are the stochastic noise 
sources / and g above which have the specified variances (11191) . To generate noises 
if>i and ip2 in Eq. [21] with variances (124H26I) we first generate two independent sets of 
random numbers with Gaussian distributions with zero mean and of variances Dq ± Do. 
Appropriate linear combinations of these random numbers yield ^i and if>2- In higher 
dimensions, stochastic forces can be analogously generated. 
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